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EVALUATION 


The  research  described  here  deals  with  the  issues  of  (1)  the  effect  on  black  box 
identification  accuracy  of  quantized  (noisy)  input/output  data,  and  (2)  the  modeling 
of  wideband  systems  by  frequency  partitioning  and  the  use  of  multirate  sampling 
within  the  sub  bands.  The  development  presented  here  shows  that  the  pencil-of- 
function  method  together  with  selected  statistical  corrections  on  the  contaminated 
data  and/or  the  use  of  multirate  sampling  leads  to  enhanced  transfer  function 
identification.  The  enhancement  is  quantitatively  described  in  terms  of  normalized 
mean  square  errors  between  the  "true"  transfer  function,  the  identified  transfer 
function  without  statistical  corrections  and  the  identified  transfer  function  with 
the  statistical  corrections. 
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I 

ADVANCED  TECHNIQUES  FOR  BLACK-BOX  MODELING 

1.  INTRODUCTION 

The  pencil-of-functions  method  is  a  black-box  modeling  method  [l]-[2]. 

Given  an  input,  output  response  pair  of  a  system  under  test,  the  algorithm 
leads  to  a  comprehensive  description  of  the  system  in  the  form  of  a  transfer 
function.  Although  the  method  was  originally  developed  for  use  upon  linear 
networks,  its  applicability  has  been  extended  by  Weiner  and  Ewen  [3]-[4]  to 
nonlinear  Volterra  models.  The  method  has  been  implemented  in  a  FORTRAN 
program  and  is  available  from  RADC  together  with  necessary  user  instructions [ 5 ] . 

The  research  described  here  deals  with  the  important  issues  of  signal  quanti- 
zation  during  analog-to-digital  conversion,  and  the  black-box  modeling  of 
wideband  systems. 

1.  Quantization 

Practical  analog-to-digital  (A/D)  converters  employ  small  word  lengths, 
typically  8  to  16  bits,  and, as  a  rule, one  can  trade  word  length  for  higher 
conversion  speed,  cost  remaining  fixed.  Unfortunately,  small  word  lengths 
lead  to  degradation  in  the  accuracy  of  the  identified  transfer  function  [6]. 

It  is  shown  here  that  the  statistical  properties  of  the  quantization  error 
can  be  exploited  to  improve  the  accuracy  ami  reliability  of  the  identified 
parameters.  The  study  thus  demonstrates  that  higher  speed  implementations 
and/or  additional  cost  benefits  may  be  achieved  for  the  pencil-of-functions 
method  than  have  heretofore  been  realized. 

2 .  Wideband  Identification 

Communication  systems  utilize  many  wideband  circuits,  for  example  , 
amplifiers  for  spread-spectrum  signals.  Black-box  modeling,  or  identification, 
of  these  circuits  poses  both  a  theoretical  and  a  practical  challenge.  A  mul¬ 
tirate  sampling  approach  to  identification  of  wideband  transmittances  is 
discussed.  It  permits  determination  of  the  transfer  function  of  a  four-to-five 
decade  bandwidth  system  from  simple  transient  tests.  Clever  selection  of 
sampling  rates  and  test  inputs  reduces  the  wideband  problem  into  three,  simpler 
smallband  problems.  The  smallband  transfer  functions  are  identified  via  the 
pencil-of-functions  method  and  then  adjoined,  systematically,  to  construct  the 
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wideband  transfer  function  estimate. 


The  report  is  structured  as  follows.  Section  2  describes  the  pencil-of- 
functions  method  in  brief.  Theoretical  details  are  omitted,  for  they  can  be 
found  elsewhere  [1],  [2].  The  description  is  included  here  for  convenience 
of  the  reader,  and  also  to  emphasize  the  discrete-time  version  of  the  method. 
A  computer  program  for  conversion  from  s  to  z  domain  transfer  functions  is 
given  in  Appendix  A.  Section  3  presents  the  study  on  improvement  of  quanti¬ 
zation-caused  degradation,  through  a  statistical  approach.  The  key  to  this 
turns  out  to  be  the  determinant  of  the  Gram  matrix  of  the  integrated  signals. 
A  computer  program,  "GQUANT",  developed  for  the  particular  case  of  impulse 
response  modeling,  is  given  in  Appendix  B.  Section  4  discusses  the  results 
of  the  study  on  wideband  systems.  Included  are  equations  and  tables  for 
ready  selection  by  the  test  engineer  of  inputs  and  sampling  rates  for  the  LF, 
MF  and  HF  band  transient  tests.  These  pulse  inputs  have  been  selected  after 
careful  study  and  are  considered  both  effective  and  laboratory  realizable. 

A  computer  program,  "USPEC",  which  generates  the  amplitude  spectra  of  the 
recommended  pulses  is  given  in  Appendix  C. 
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2.  PENCIL-OF- FUNCTIONS  METHOD 


Recorded  input,  output  responses  of  a  network  can  be  integrated  to  yield 
a  family  of  signals,  called  measurement  signals.  Application  of  the  pencil- 
of-functions  theorem  [1]  to  this  family  yields,  in  a  closed  form,  the  identi¬ 
fied  parameters  of  the  network  function.  The  procedure  for  this  black-box 
modeling  method  is  described  below.  Although  proofs  are  omitted,  the  usefulness 
of  the  technique  will  be  demonstrated  with  examples.  Discrete-time  signals 
are  chosen  for  the  presentation  here,  because  of  inherent  computational  advan¬ 
tages,  although  such  signals  must  often  be  obtained  by  sampling  a  continuous¬ 
time  system. 

2.1  SIMULTANEOUS  NUMERATOR  AND  DENOMINATOR  DETERMINATION 
Identification  Problem 

Given  the  input-output  observations 

{ u (k) } ,  { y (k) },  k=0, 1 , . .  ,K  (1) 

arising  from  a  physical  system  (see  Fig.  1)  believed  to  be  linear,  and  of  finite 
order,  it  is  desired  to  find  a  system  model 


H  (z) 


V’1 

+  .  .  .  . 

+ 

G 

1 

-1 

-n 

1  +  a  .z 

+  .  .  .  . 

+ 

a  z 

1 

n 

n  d  .  z 


l 


(2) 


(3) 


which  best  fits  the  observations,  in  some  sense  (see  Fig.  2).  A  solution  can 
be  obtained  by  use  of  the  pencil-of-f unctions  theorem  as  discussed  below. 

For  convenience  denote  sequences  (u(k)}  and  { y (k) }  simply  as  u  and  y, 
respectively.  Also,  denote  the  inner-product  of  two  sequences  as 

def  K 

<  x ,y  >  =  E  x(k)  y(k)  (A) 

k=0 


Measurement  Sequences 

Krcm  the  given  sequences  y  and  u  we  form  the  following  set  of  sequences, 
rilled  measurement  sequences: 


3 


_ ,  1  ^ 


Donverte 
S.I.=  A 


| _ |  u(kA)=u(k) 

Fig.  1.  Response-pair  from  system  under  test 


Fig.  2.  Identification  problem:  Find  H(z)  so  that 
y(k)  is  close  to  y(k) 
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y-^k)  =  y(k) 

y,(k)  =  y1(0)+y1(l)+. . .+y1(k) 

•  •  •  •  •  •  • 

a  a  a  a  a  a  a 

yn+i(k)  =  yn(0)+yna)+--*+yn(k)  (5) 


ux(k)  =  u(k-l) 

u2(k)  =  u^OJ+u^d)  ....  +u^ (k) 

u  (k)  =  u  (0)+u  (1)+....  +u  (k)  (6) 

n+1  n  n  n 

where  n  is  the  order  of  the  model  desired.  That  is,  n  is  the  degree  of  the 
network  function  H. 

Note  that  these  sequences  represent  repeated  discrete  integrations  of 
the  observed  signals  y(k)  and  u(k),  respectively,  i.e., 


Vi(k) 

k 

=  E  y.OO 

1=0  3 

j=l . 

n 

(7) 

Vi00 

k 

=  E  u. (£) 

«,=0  3 

j=l»  •••» 

n 

18) 

Iquivnlently,  yj+^(k)  is  obtained  by 
,/(z-l)  as  shown  in  Fig.  3.  Likewise 

passing  y  (k) 
,  u  .j  (k)  is 

through 

obtained 

the  filter  I(z) 
by  passing  u^ (k) 

through  the  discrete  integrator  l(z). 

Gram  Matrix 

Next  form  the  following  inner-product  matrix 


5 


Digital  integrator  I(z)= 


z  -  1 


Fig.  3.  Generation  of  integrated  signals  and  the  Gram  matrix 
(u(k)  included  in  the  inner-product  generator  if 
model  is  required  to  have  direct  transmission  ) 


where  we  have  used  the  notation  N  =  n+1  for  convenience.  This  (N+n)  x  (N+n) 
dimensional  matrix  is  the  Gram  matrix  [10]  of  the  (N+n)  dimensional  vector 
sequence 


{f (k)}. 


k  =  0,1, ,  K 


where 


f(k)  = 


yL(k) 

y2<k) 


yN(k) 
,  u2(k) 


UN^> 


To  state  this  observation  formally,  we  have 

K  T 

F  =  Z  f (k)f  (k) 

k=0 


(10) 


(11) 


(12) 


^The  entrv  u^(k)  is  omitted  in  _f(k),  and  therefore  in  the  formation  of  the 

gram  matrix  F,  whenever  direct  transmission  in  the  model  is  absent  (that  is 
when  the  coefficient  b  in  the  function  H(z)  is  constrained  to  be  zero). 


r  t 


Diagonal  Cofactors 

Denote  the  diagonal  cofactors  of  F  as  D^: 

D^  =  i,i  cofactor  of  F  (13) 

Recall  that  the  i,i  cofactor  of  a  square  matrix  is  the  determinant  of  the 
matrix  after  deleting  the  ith  row  and  the  ith  column. 

Parameters  of  the  Network  Function 

The  parameters  of  the  network  function  are  given  by  the  square-roots 
of  D^  up  to  a  multiplicative  constant.  That  is 

N  -  i-1  n  1  i-1 

[  E  /D.  (1  -  z_i)  ]  Y (z)  =  [  E  /D  .  z_i(l  -  z'1)  ]  U (z)  (14) 

.  .1  .  ,  N+x 

1=1  i=l 

which  can  be  normalized,  by  dividing  by  D  =  /D  +  . . .+*4)^,  so  that  the  leading 

coefficient  becomes  unity.  Clearly  the  computed  transfer  function  becomes 

n  ,  i-1 

2  (  E  /DN+.(l-z  X)  ]  /D 

H  (z)  = - — -  (15) 

N  _  i-1 

[  £  t/0.  (1-z  x)  ]  /D 

i=l  1 


REMARKS 

•  Note  that  the  first  measurement  signal  is  the  network  output  itself,  y^  =  y. 
Next  follow  its  successive  integrations.  Each  of  these  signals  can  be 
expressed  directly  in  terms  of  y(k).  Indeed,  if  we  let  I(z)  =  z/ (z-1) ,  Y^+^(z) 
=  1^ (z)  Y(z)  so  that  Yj+^(k)  =  i ^ (k)  ^*)  y(k)  where  i  (k)  is  the  inverse 
transform  of  1^ (z)  and  ©  denotes  discrete-time  convolution. 

•  The  dimensionality  of  the  measurement  vector  f(k)  is  2n+l  =  N+n  when  the 

direct  transmission  term  b  is  constrained  to  be  zero.  If  the  network  does 

o 

have  direct  transmission,  u^(k)  =  u(k)  should  be  included  in  the  vector  f_ 
so  that  its  dimensionality,  as  well  as  that  of  the  corresponding  Gram  matrix 
F,  becomes  2n+2  =  2N.  The  right  hand  side  of  equation  (14)  modifies  slightly 
as  follows 

N  i-1  N  i-1 

[  Z  /D.  (1  -  z_i)  ]  Y(z)  =  [  E  /D  (1  -  z  X)  ]  U(z)  (16) 

i-1  1  i-1  N+1 

The  counterpart  of  equation  (15)  follows  from  (16)  and  is  therefore  not  given 
here. 


_ •  s 
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To  illustrate  the  steps  of  the  method,  a  simple  example  is  given 
next.  (The  reader,  unfamiliar  with  the  pencil-of-f unctions  method,  may  wish 
to  work  the  details  with  pencil  and  paper;  others  may  skip  this  example.) 


Examp le  1 

Consider  the  setup  of  Fig.  4  where  u^(k)  denotes  the  input  signal  and 
y^(k)  the  output.  The  network  is  known  to  have  direct  transmission  and  of 
first  order  (i.e.,  the  s-domain  transfer  function  is  of  the  type  (d^s  +  d^)/ 
(s  +  Cq)).  The  measurements  are  made  every  1  ms  for  5  samples,  k  =  0,1,..., 4. 
Unit  pulse  input 

Suppose  the  following  signals  are  generated  as  a  result  of  a  unit  pulse 


(only  y^  and  u^ 

may  have  been  recorded  in  real 

time) : 

Vj  (k) 

1.0 

1.2  0.96 

0.768 

0.6144 

y,00 

1.0 

2.2  3.16 

3.928 

4.5424 

ux(k) 

1 

0  0 

0 

0 

u2(k) 

1 

1  1 

1 

1 

The  Gram  matrix 

of  the  signals  y^,  y0 

,  u^  and  u2 

is 

F  = 

_ 4. 3289 

12.4811 

1.0 

4.5424 

12.4811  1.0 

51.8881  1.0 

1.0  1.0 

14.8304  1.0 

4.5424 

14.8304 

1.0 

5.0 

which  yields  the  following  square-roots  of  the  diagonal  cofactors. 

/D,  =  3.5032  /D„  =  0.87581  /D_  =  1.7516  /D.  =  -6.1307 

12  3  4 

Note  that  the  signs  of  these  square-roots  are  chosen  in  direct  correspondence 

with  the  signs  of  the  cofactors  of  the  first  row  of  F  [1].  Now,  substitution 

into  (16)  and  division  by  (D^  +  D^)  leads  to  the  equation 

(1  -  0.8  z_1)  Y(z)  =  (1  +  0.4  z"1)  U (z) 

Clearly,  the  true  parameters  have  been  recovered. 
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z,  .  _ 2.00z~2  -  3.71150Z  3  +  1.7128z  4 _ 

H  -1  -2  -3  -4 

1  -  3.3792z  +  4.4286z  -  2.7181z  +  0.66898z 

s-poles :  (-0.002  +  j  0.0699714) (106) 

(-0.399  +  j  1.131373  )O06) 

Using  the  inverse  of  the  pole- zero  transform,  the  s-domain  transfer 
function  can  be  obtained.  The  poles  turn  out  as  shown  above. 

The  response  of  the  model  and  the  actual  network  response  are  compared 
in  Fig.  5b. 
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microsec. 


REMARKS 


•  When  the  network  under  test  is  of  order  n,  i.e.,  when  the  model  order 
is  equal  to  the  intrinsic  order  of  the  network,  the  rank  of  the  matrix  F 
equals  its  dimensionality  minus  one. 

•  The  matrix  F  is  positive  semi -definite. 

•  In  actual  application  the  matrix  F  will  be  formed  from  quantized  versions 

of  signals  y  and  u.  Call  this  corrupted  matrix  as  G.  It  will  be  shown  that 
2 

E{g}=F  +  a  P,  where  P  denotes  the  correlation  matrix  of  unit  noise  and  E  denotes 
the  statistical  expectation  operator.  It  will  be  shown  in  Section  III  that  E{g} 
has  full  rank  (equal  to  the  dimensionality  of  F) . 

As  seen  earlier,  the  pencil  of  functions  method  uses  the  square-roots  of 
the  diagonal  cofactors  of  F.  A  very  important  advantage  of  the  method  is  the 
following. 

"Since  F  is  positive  semi-definite  (G  positive  definite  with  unit  prob¬ 
ability),  its  diagonal  cofactors  are  non-negative  (strictly  positive).  Hence, 
there  is  a  built-in  check  and  stopping  point  when,  due  to  computational  errors 
or  wrong  choice  of  model  order,  one  or  more  of  these  cofactors  turns  out  to 
be  negative." 

The  computations  involve  finding  the  cofactors  of  a  2n  +  1  or  2n  +  2 
dimensional  matrix.  For  the  special  case  of  impulse  response  modeling  the 
calculation  of  denominator  and  numerator  coefficients  can  be  decoupled,  so 
that  computations  involve  only  an  n  +  1  dimensional  matrix.  This  will  be 
discussed  next  . 


2.2  DECOUPLED  PROCEDURE  FOR  MODELING  IMPULSE  RESPONSES 
Consider  that  y(k)  is  the  impulse  response  of  a  network  and  that  a 
suitable  K  has  been  selected  such  that  y(k)  -  0  for  k>K.  We  define  the 
reverse-time  integrated  signals  as  follows  [2],  [11] 


yL(k)  =  y (k) 

y2(k)  =  y1(k)+. . . .+y1 (K) , 

yN(k)  =  yn(k)+. . . .+yn(K), 


(17) 
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(Recall,  N=n+1) .  Let  F  be  defined  as 


<y 


1’ 


v  • 


F 


V 


<yl* 


<yN’V 


^i’v 


K 

E  y.(k)y.(k) 
k=l  1  1 


(18) 


or,  equivalently, 

K  T 

F  =  E  f_(k)f. (k) 

k=l 


(19) 


where  f_T(k)  =  [y1(k)  y2(k)  -  yN(k)]. 


Then,  it  can  be  shown  that  the 


denominator  polynomial  is  given  by 
N  . 

A(z)  =  z~n[  E  /D.  (z-1)N_1]//D 
i=l  1 


(20) 


where  denotes  the  ith  diagonal  cofactor  of  the  matrix  F.  Note  the 
positive  powers  of  z  on  the  right  hand  side.  Further,  the  numerator  co¬ 
efficients  are  obtained  as: 


f'b  ~] 

°  ! 

P11  P1N 

P  V 

t 

. 

I 

■  !  _ 

•  *  *  *  *  * 

1 

•  1 

. 

1 

b  j 

L  n  J 

_  PN1 -  PNN  _ 

L  qN_ 

(21a) 


p^.  =  <w(k+L-i),  w(ktl-j)> 

q..  =  <  y (k) ,  w(k+l-i) > 

where  w(k)  is  the  impulse  response  (i.e.,  inverse  z-transform) 
and  inner  products  are  summed  from  k=0  to  K. 


(21b) 

(21c) 
of  1/A(z) 


If  the  network  is  known  to  have  no  direct  transmission,  i.e.,  j.s  suspected 
to  be  zero,  then  N  should  be  replaced  by  n  on  the  right  hand  side,  b  by 
b  and  in  forming  the  inner  products  w(k+l-i)  should  be  replaced  by  w(k-i) 
(likewise,  w(k+l-j)  should  be  replaced  by  w(k-j)). 

Three  examples  will  be  presented  next.  The  first  is  a  simple,  paper- 
pencil  type  example;  it  considers  the  same  impulse  response  as  did  Example  1 
(page  9)  but  with  a  long  record  length.  The  final  example  is  interesting 
because  it  deals  with  an  impulse  response  which,  theoretically,  requires 
an  infinite  order  system  (of  type  1))  for  exact  reproduction;  a  fifth 
order  model  is  computed  by  the  pencil-of-f unctions  method  which  yields  a 
fractional  energy  error  of  0.0359.  In  the  first  two  examples  the  true 
transfer  function  is  recovered  by  the  modeling  technique,  i.e.,  the  frac¬ 
tional  energy  error  is  zero. 

Notation  - 

y (kA)  or  y (k) 

y  (k  A)  or  y  (k) 

K  2 

S  =  E  y  (k) 
k=  0 

^  ^  2 

e  =  Z  y  (k) 
k=0 

v  =  e  /S 
n  =  loo(i-v) 

Example  3 

Given  the  left  hand  side  of  y^(k)  =  1.5(0. 8)^-  0.56^,  we  find  for  K=40 


Model  response 

Model  response  error  y(k)-y(k) 

Response  energy 

Error  energy 

Fractional  energy  error,  or  simply  fractional 
error,  or  normalized  mean  square  error 

Per  cent  modeling  efficiency 
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!  4  20’1 

_  _  j  ;  (inner  products  are 

I  from  k=l  to  40) 

_20  100J 

Then  =  100  and  =  4.  Equation  (20)  yields 

A(z)  =  z_1 (10z-8) /10  =  1  -  0. 8z_1 

Equation  (21),  in  turn  becomes 


which  produces  B(z)  =  (1  +  0.4z  ).  The  model  has  been  identified  perfectly 

with  zero  fractional  error. 

Example  4 

A  fourth  order  network  is  known  to  have  zero  direct  transmission  (bQ=0) . 
The  numerical  data  of  its  impulse  response, 

y(t)  =  10  e~2t  Sin (2t)  -  2  e"°‘5t  Sin(4t) 
is  recorded  at  uniformly  sampled  intervals  of  A  =  0.2  sec.  For  K=150  (which 
signifies  a  long  record;  K  =  30  sec),  we  find 


4.51403  2.38856 

0.779798 

0.114119 

-0.05724740 

2.09775 

0.953032 

0.253410 

0.00074006 

F  = 

det 

F  =  0.54E-14 

0.501610 

0. 157312 
0.057232 

0.01499790 

0.00960080 

0.00352616 

Note  -  All  summations  have  employed  a  multiplication  factor  A,  for  scaling 
purposes,  both  in  forming  the  integrated  signals  and  in  forming  the  inner- 
produets.  However,  to  undo  the  effect  of  this  scaling,  the  ith  diagonal 
cofartnr  has  to  be  multiplied  with  A2^  to  yield  D..  The  entire  process  will 
be  called  .'.-scaling 

The  values  of  vD./D.  are 
1  1 

1  1.50410  1.33762  0.58517  0.11959 
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Equations  (20)  and  (21)  yield  the  following  denominator  and  numerator 
coefficients. 

Denominator  1  -2.49588  2.82521  -1.5760  0.36786 

Numerator  0  -1.31238  1.68950  -1.55568  0.00158 

The  fractional  energy  error  turns  out  to  be  V=  0.1E-6.  As  seen  from 

A 

Fig.  6  the  model  response  y  is  indistinguishable  from  the  true  response  y. 
Example  5 

Here  we  consider  a  problem  in  approximation.  This  terminology, 
rather  than  identification,  is  appropriate  since  the  square  pulse 

y(t)  =  f  1  for  0  <  t  _<  1 

£  0  otherwise 

cannot  be  exactly  reproduced  as  the  impulse  response  of  a  finite  order 
linear  system.  A  fifth  order  model  is  desired  whose  impulse  response 
approximates  this  signal.  Using  A=0.05  sec.  and  setting  y(0)=0,  y(k)=l 
for  k=l,..,20,  y(k)=0  for  k=21,..,40  the  following  Gram  matrix  is  obtained. 


F 


1.0  0.525000 

0.192500 

0.0553438 

0.01328250 

0.002767190 

0.358750 

0.146781 

0.0635731 

0.0448284 

0.0201175 

0.0065189 

0.01117940 

0.00514064 

0.00169444 

0.00044604 

0.002391640 

0.001119480 

0.000373780 

0.000099363 

0.000022306 

det  F  =  0.344E-26  (note  -  A  scaling  is  employed) 


The  values  of  A)./D,  are 

l  1 

1  1.1458  0.67865  0.22900  0.042415  0.003383 

Equations  (20)  and  (21)  yield  the  following  z  transfer  function  coeffi¬ 
cients  . 

Denom.  1  -3.854184  6.095388  -4.93206  2.037106  -0.0342867 

Numer.  0  -1.179760  3.806755  -5.19377  3.53997  -1.039771 

The  fractional  energy  error  turns  out  to  be  V=0,0359  with  a  corresponding 

A 

modeling  efficiency  of  96.4%.  The  model  response  y  is  compared  with 


the  desired,  ideal,  response  y  in  Fig.  7a. 

A  brute  force  application  of  the  correction  procedure  given  in 
the  next  section  (Section  3)  results  in  the  model  response  shown 
in  Fig.  7b.  Of  course  better  approximations  can  only  be  obtained  with 
higher  order  models. 
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3.  QUANTIZATION  ERROR:  IMPROVEMENT  OF  ESTIMATES 
(PENCIL-OF- FUNCTIONS  METHOD) 


Practical  analog-to-digital  (A/D)  converters  employ  small  word  lengths, 
typically  8  to  16,  and  as  a  result  incur  quantization  error  in  the  repre¬ 
sentation  of  the  signal.  This,  in  turn,  causes  degradation  in  the  accuracy 
of  the  identified  transfer  function  [6].  It  will  be  shown  in  this  section 
that  the  statistical  properties  of  the  quantization  error  can  be  exploited 
to  improve  the  accuracy  of  the  parameter  estimates.  A  computer  program 
"GQUANT"  incorporating  the  technique  developed  is  given  in  Appendix  B. 

The  principle  of  analog-to-digital  conversion  is  explained  well  in 
references  [6],  [9].  For  our  purposes  certain  essential  properties  are  most 
pertinent.  If  b  bits  are  used  (including  the  sign  bit)  and  XMSB  is  the  analog 
value  of  the  most  significant  bit  (next  to  the  sign  bit),  then  the  following 
observations  and  properties  hold. 

(a)  The  step  size  equals 


6  = 


XMSB 

nb-2 


(22) 


(b)  For  an  input  y  to  the  A/D  convertor  the  analog  value  of  the  output  is 


x  =  y  +  e 

6 

where  |e|<  for  roundoff  and  |e|<  5  for  truncation  . 


(23) 


(c)  If  the  signal  excursions  during  each  sampling  time- interval  A 
are  large  compared  to  6,  then 

x(k)  =  y (k)  +  e (k)  y(k)d2f  y(kA)  (24) 


where  e(k)  is  an  independent  sequence  of  random  variables  having 
a  uniform  distribution  over  one  step  size  6.  In  case  of  roundoff  , 
this  distribution  is  centered  at  zero,  so  that  the  random  variable 
e(k)has  a  zero  mean  and  a  variance [7]. 

,2 

Var{e(k)  }  =  (25) 

In  the  ensuing  discussion  we  will  assume  the  A/D  converter  employs 
roundof  f . 

(d)  Under  the  assumptions  in  (c)  above,  the  error  sequence  e(k)  is 
uncorrelated  with  the  parent  sequence  y(k) 
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Simulation  shows  that  neither  of  the  properties  (c)  or  (d)  strictly 
hold  in  practice.  However,  we  will  use  these  properties  exercising  caution 
where  necessary. 

For  definiteness  we  will  discuss  in  detail  the  correction  technique  for 
impulse  response  modeling  method  of  subsection  2.2.  Parallel  formulas  are 
applicable  to  the  simultaneous  denominator  and  numerator  modeling  procedure 
of  subsection  2.1,  but  will  not  be  given  here.  Recall  that  the  poles  of  the 
model  are  obtained  from  the  Gram  matrix  of  the  signal  y  and  its  successive 
integrations.  We  therefore  begin  with  the  analysis  and  correction  of  the 
quantized  Gram  matrix. 

3.1  GRAM  MATRIX  OF  THE  QUANTIZED  SIGNAL 

We  will  use  the  model  of  equation  (24)  for  the  quantized  signal 
x(k)  =  y(k)  +  e(k) 

where{E  e(k)}=  0,  E{e(k)  e(£)}=  0  and  E{y(k)  e(£)}=  0  for  all  k  and  £.  For 
the  reversed  time  integrated  signals  we  have 

x^k)  =  yt(k)  +  e1(k)  (26) 

where  e^(k)  are  derived  from  e(k)  through  equations  analogous  to  (17),  i.e.. 


ei+1(k)  =  ei(k)  +  e^k+1)  + 
Define  also  the  vector  sequences 

PX1 (k)l 

!  x2(k)  | 


+  e±(K) 


x  (k) 
.  N  . 


£.00  = 


e1(k) 

e,(k) 


e  (k) 
.  N  i 


,  k=0, . . . , K  (28) 


Then  the  Gram  matrix  of  the  quantized  signal  can  be  written  as 


G  =  E  £(k)  gT(k) 


=  Z  [f(k)  fT(k)  +  f(k)  pT(k)  +  P(k)  fT(k)  +  P(k)£T(k)] 


Observation  1 


E  (G>  =  E  f(k)fT(k)  +  E{  E  £(k)£T(k)} 

k=0  k=0 

=  F  +  o2  P  (30) 

where  P  is  the  unit  noise  covariance  matrix  defined  below.  Further,  if 
properties  (c)  and  (d)  strictly  hold  then 


Observation  2 

The  unit  noise  covariance  matrix  is  given  by 
K  T 

P  =  E  {  E  ja_(k)  P  (k))  (31) 

k=0 

where  P(k)  =  [e^(k)  e^(k)  ..  e^(k)]^  as  before,  but  ej  (k)  =  e(k)  is  taken 
to  be  a  zero  mean,  unit  variance,  uncorrelated  sequence. 

Remark 

2 

If  properties  (c)  and  (d)  do  not  strictly  hold,  then  the  value  of  O 
(and  possibly  the  definition  of  P)  should  be  modified.  We  will  estimate  o2 
so  as  to  satisfy  the  following  criterion. 

Jain's  Identification  Criterion 

Consistent  with  the  noise  and  bias  models  the  estimated  Gram  matrix 
should  achieve  a  minimum  of  the  determinant . 

Whatever  method  is  used  to  choose  the  estimated  Gram  matrix,  care 
should  be  taken  to  make  sure  that  its  determinant  remains  nonnegalive, 
since  the  determinant  of  the  true  Gram  matrix  is  nonnegative  (see  page  13  ). 
An  approach  to  estimation  of  the  Gram  matrix  is  presented  in  subsection  3.4. 
First,  however,  we  discuss  the  computation  of  the  unit  noise  covariance 
matrix . 

3.2  UNIT  NOISE  COVARIANCE  MATRIX 
Examination  of  the  sequences 

Oj (k)  =  e (k) 

e2(k)  =  e (k)+  . . .  +  e(K) 

e3(k)  =  e(k)  +  ...  +  (K+l-k)e(K) 

e4(k)  =  e(k)+  ...  +  (K+l-k) 2e (K) 
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leads  to  the  general  formula 


IN.  .  - 

e...^)  =  Z  £  e(k-£)  (33) 

1  £= 0 

We  then  have  (using  the  definitions  in  (28)  and  (31)) 

K  1 

!  i  1 

K-k  ' 

p(k)  =  Z  .  £  ;  e(k-£)  (34a) 

1=0 


-£n"lJ 

K-k 

=  Z  r(£)  e (k-£)  (34b) 

£=0 

where  r*(£)  =  [6^  1  ...  £n_i];  <5^  is  the  unit  pulse  sequence.  Then 

K  K— k  K-k 

P  =  E{  Z  Z  Z  r_(H)  _rT  (m)  e(k-£)  e(k-m)} 

k=0  £=0  m=0 

and,  since  e  is  a  zero  mean,  unit  variance,  uncorrelated  sequence, 

K  K-k  T 

P  =  Z  Z  r(£)  _r  (£) 
k=0  £=0 

=  Z  (K-k+1)  r(k)  rT(k)  (33) 

k=0 

Note  that  P  is  determined  entirely  by  the  integers  N  and  K,  the  dimensionality 
of  F  (recall  N=n+1)  and  the  length  of  the  observed  sequences,  respectively. 
Clearly,  P  can  be  precomputed  and  stored. 


3.3  ESTIMATION  OF  QUANTIZATION  ERROR  VARIANCE 
The  discussion  in  subsection  3.1,  specifically  equation  (30),  leads  us 
to  estimate  F  as 

F  =  G  -  02P  (36) 


TT 
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where  0  will  be  chosen  so  as  to  minimize  the  determinant  of  F. 

One  possible  approach  to  this  minimization  is  developed  here.  We  use 
the  fact  that  the  rank  of  the  true  Gram  matrix  F  is  n,  i.e.,  its  determinent 
is  zero.  Rewriting  (36) 

-  2 

F  =  G  -  a  P 


we  set  the  determinant  of  both  sides  to  zero.  If  the  quantization  error  is 
small,  we  can  approximate  the  determinant  of  the  right  hand  side  by  the  first 
two  terms  of  the  determinant  expansion  theorem.  Thus 


| F |  =  0  -  G-O2  £  det[G,P]  (37) 

where  the  notation  [G,P]^  means  the  matrix  obtained  by  replacing  the  ith 
column  of  G  by  the  ith  column  of  P. 

Then 


*2  _  I G 1 

Sdet[G,P], 


(38) 


and,  of  course, 

F  =  G  -  OP 


(39) 


Note  that  formula  (39)  can  also  be  applied  recursively,  by  replacing  G  in  (39) 
with  the  last  estimate  of  F.  An  exit  must  be  made  when  the  determinant  of  the 
estimated  matrix  ceases  to  reduce  further  (or  begins  to  increase). 


3.4  SIMULATION  EXAMPLES 

As  stated  earlier,  a  FORTRAN  IV  computer  program  "GQUANT"  has  been  developed 
for  simulation  and  modeling  of  quantized  impulse  responses.  A  rational  mo.ul  of 
the  type  given  in  equation  (1)  is  produced,  except  that  b^  is  constrained  to 
zero;  i.e.,  the  network  is  assumed  to  have  no  direct  transmission.  (Slight 
modification  in  the  computation  of  numerator  coefficients  enables  this  con¬ 
straint  to  be  removed.)  Equivalent  s-domain  description  can  be  obtained  through 
appropriate  z  to  s  transformation.  Salient  features  of  the  program  are  the 
following. 

It  can  be  used  in  either  a  simulation  mode  (IRESP=1  or  2)  or  in  network- 
response-data  entry  mode  (IRESP=0) 

Model  can  be  obtained  for  unquantized  signal  (ISPN=  -1,  IFIX  =-l» 

NFTX  immaterial)  when  in  the  simulation  mode,  or  actual  response-data  when 
in  data  entry  mode. 
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Model  can  be  obtained  for  the  quantized  signal  (ISPN“1)  without  any  sta¬ 
tistical  correction  (IFIX=  -1).  Intended  for  use  in  simulation  mode. 
Model  can  be  obtained  for  the  quantized  signal  (ISIM=1  or  2  and  ISPN=1, 
or  ISIM=0  and  ISPN=0)  with  statistical  correction  (IFIX=1) ;  the  use  of 
IBIAS=1  performs  a  bias  correction  in  addition  to  statistical  noise 
correction. 

Two  examples  are  given  below,  one  in  which  a  second  order  network  response 
is  simulated  and  another  in  which  a  fourth  order  response  is  simulated.  Thus 
both  examples  pertain  to  simulated  impulse  responses. 


H.  ,  , (z)  True  transfer  function  of  the  network* 

ideal 

H(z)  Transfer  function  obtained  by  application  of  pencil-of- 

functions  method  upon  unquantized  signal.  Note  that  H(z) 

need  not  be  equal  to  H.,  (z);  among  the  reasons  for 

ideal 

this  are  computation  errors,  and  the  use  of  K  ^ 

Hquant  Transfer  function  obtained  by  from  the  quantized  signal. 

(no  correction  is  applied) 

H(z)  Transfer  function  obtained  from  the  quantized  signal; 

one  or  more  iterations  of  statistical  correction  for 
^  quantization  errors  are  used. 

H(z)  Transfer  function  obtained  from  the  quantized  signal;  in 

addition  to  statistical  correction  for  quantization  errors, 
correction  is  also  applied  for  possible  bias  in  the  data. 

It  should  be  mentioned  that  the  usefulness  of  bias 
correction  arises  both  because  the  quantization  errors  in 
particular  record  of  data  may  not  be  zero-mean,  and  also 
because  K  4  00  may  produce  an  apparent  bias  in  data. 

NDIG  Length  of  binary  word  b^^  ...  (note  bNDIG  is  the 

sign  bit,  bNDIG_j_  the  most  significant  bit,  ...,  and  b^ 
the  least  significant  bit;  also,  we  have  employed  a 
mid-tread  type  of  quantizer  in  simulation) 

The  analog  weight  (or  value)  of  the  most  significant 
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Example  6 

A  second  order  network  with  zero  direct  transmission  is  simulated.  Its 
impulse  response 


,.s  0  ~2t  *  -0 . 5t 

y (t)  =  2  e  -  2  e 

is  sampled  uniformly  at  intervals  A  =  0.2  sec.  apart. 

transfer  function  H.,  , (z)  are 

ideal 


The  coefficients  of  the 


Denominator 

1  -1.575157  0.606530 

Numerator 

0  -0.469035  0 

Without  quantization  the  modeling  program  yields  (using  ISPN  =  -1, 
IFIX  =  -1,  NFIX  immaterial)  the  following  results: 

H(z)  (using  ISPN  =  -1,  IFIX  =  -1,  NFIX  immaterial)  V  =  0.6E-8 
Denominator 

1  -1.575180  0.606551 

Numera c  or 

0  0.469050  -0.000041 

Experiment  1 

For  XMSB  =  5.0  Volts  and  NDIG  =  10,  the  program  yields  the  following 
results 


H  (z)  (using  ISPN  =  1,  IFIX  =  -1,  NFIX  immaterial)  V  =  0.84E-3 

quant 

Denominator 

1  -1.635921  0.662757 

Numerator 

0  0.513516  -0.114095 

H(z)  (using  ISPN  =  1,  IFIX  =  1,  NFIX  =  3;  includes  bias  correction)V=0.62h-3 
Denominator 

1  -1.628129  0.655569 

Numerator 

0  0.506869  -0.098299 

The  impulse  responses  of  ^  ant(z)  and  H(z)  are  compared  with  that  of 

^ideal^2^  '*'n  ^8*  ®a  nnc*  (The  quantized  signal  used  in  determining  these 

transfer  functions  is  shown  in  Fig.  8c.)  Although  the  improvement  through  sta¬ 
tistical  correction  is  hard  to  discern  from  these  figures,  the  fractional  energy 
i  rrur  clearly  points  to  a  slight  improvement. 
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Fig.  8a  Comparison  of  network  and  model  impulse  responses. 
Mode)  is  obtained  from  quantized  data. 

(No  statistical  correction  is  applied) 

-2  8- 
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Fig.  8b  Comparison  of  network  and  model  impulse  responses. 

Model  is  obtained  from  quantized  data;  statistical 
correction  is  applied. 
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A  more  impressive  improvement  is  achieved  in  the  next  experiment. 

Experiment  2 

For  XMSB  =  5.0  Volts  and  NDIG  =  7  ,  the  program  yields  the  following 


results 

H  (z)  (using  ISPN  =  1,  IFIX  =  -1, 

quant 

NFIX  immaterial) 

V  =  0.0061 

Denominator 

1 

-1. 721641 

0.743261 

Numerator 

0 

0.571452 

-0.263324 

H(z)  (using 

ISPN  =  1,  IFIX  =  1,  NFIX 

=  1;  IBIAS  =  0) 

V  =  0.0047 

Denominator 

1 

-1. 703728 

0. 726890 

Numerator 

0 

0.550597 

-0.219990 

H(z)  (using 

ISPN  =  1,  IFIX  =  1,  NFIX 

=  1;  IB  LAS  =  1) 

v  =  0.0042 

Denominator 

1 

-1.696989 

0. 720718 

Numerator 

0 

0.543307 

-0.20437 L 

Example  7 

A  fourth  order  network  with  zero  direct  transmission  is  simulated.  the 
impulse  response 

y(t)  =  10  e  Sin(2t)  -  2  e  sin(4t) 

is  sampled  uniformly  at  intervals  A  =  0.2  sec.  apart.  The  coefficients  of 

the  transfer  function  H.  ,  ,(z)  are 

ideal 


Denominator 

1  -2.495629  2.824925  -1.577498  0.367879 


Numerator 

0  1.312168  -1. 688] 52  1.553863  0 

Without  quantization  the  modeling  program  yields  (using  1SPN  =  -1, 
IFIX  -  - 1 ,  NFIX  immaterial)  the  following  results: 
iifz)  tuning  ISPN  =  -1,  IFIX  =  -1,  NFIX  immaterial)  v  =  0.1F.-6 


3] 


Denominator 

1  -2.495883  2.825209  -1.577598  0.367863 

Numerator 

0  -1.312376  1.689499  -1.555676  0.001578 

For  XMSB  =  5.0  Volts  and  NDIG  =  10,  the  program  yields  the  following 
results 


Hquant(z)  (usin8  ISPN  =  IFIX  =  -1.  NFIX  immaterial)  V  =  0.076 


Denominator 

1 

-3.003323 

3.619004 

-2.057970 

0.455114 

Numerator 

0 

-1.290046 

2.315347 

-1.991812 

0.834767 

H(z)  (using 

ISPN  =  1,  IFIX 

=  1,  NFIX  = 

1;  IBIAS  =  0) 

V  =  0.045 

Denominator 

1 

-2.930222 

3.483747 

-1.979050 

0.440495 

Numerator 

0 

-1.368196 

2.599335 

-2.567285 

1.185892 

H(z)  (using 

ISPN  =  1,  IFIX 

=  1,  NFIX  = 

l;  IBIAS  =  1) 

V  =  0.040 

Denominator 

1 

-2.910943 

3.447190 

-1.956858 

0.436115 

Numerator 

0 

-1.393888 

2.690009 

-2.731991 

1.280967 

Clearly, 

a  reduction  in 

energy  error  has  been  achieved  via  statistical 

correction. 

Remarks 

The  application  of  the 

statistical 

correction  was 

predicated  upon  several 

assumptions.  Experiments  show  that  these  assumptions  are  not  satisfactorily 
met.  The  following  comments  therefore  arise. 

The  quantization  error  process  e(k)  is  not  white.  It  might  be  useful 
in  future  work  to  model  this  error  process  as  a  first  order  process  and  estimate 
the  corner  frequency  of  this  process  together  with  its  intensity. 

•  The  correlation  between  the  quantization  error  e(k)  and  the  input 

signal  y(k)  is  not  zero.  This  may  be  ameliorated  by  the  use  of  a  well  known 
technique  [13]  namely  the  addition  and,  after  quantization,  the  subtraction 
of  a  dither  signal^.  This  is  shown  in  Fig.  9  .  The  application  of  this 

2Pseudo-random  binary  signals  are  often  used  as  dither  signals. 
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technique  to  our  problem,  and  the  extent  of  improvement  achieved [ 14 ] ,  remain 
subjects  of  future  investigation. 


Known  dither 
signal 

Fig.  9  Use  of  dither  signal  to  decorrelate  y(k)ande(k) 

•  In  estimating  the  intensity  (variance)  of  noise  via  equation  (36) 

only  the  first  two  terms  of  the  determinant  expansion  were  retained.  Perhaps 

three  terms,  i.e.,  constant  linear  and  quadratic,  should  be  retained  in  order 

2 

to  get  a  more  accurate  estimate  of  O  .  However,  we  feel  that  the  benefit  of 
this  step  would  be  realized  only  after  the  steps  1  and  2  stated  above  have 
been  taken. 
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4.  WIDEBAND  IDENTIFICATION 

Determining  the  transfer  function  of  a  network  from  its  ob¬ 
served  input-output  responses  represents  the  inverse  of  the  analysis 
problem.  Interest  in  this  problem  arises  from  the  frequent  need  for 
a  relatively  simple  mathematical  description  of  the  system  so  that 
behavior  for  other  anticipated  inputs  may  be  predicted  up  to  accept¬ 
able  accuracies.  However,  the  identification  of  wideband  networks 
presents  some  unique  difficulties.  Consider,  for  example,  a  network 
whose  frequencies  of  interest  range  from  fQ  Hz  to  (10^) f  Hz.  To 
identify  the  corner  frequencies  at  the  low  end,  one  would  require  an 
observation  record  of  length  T  -  1/f  sec.  On  the  other  hand, 
in  order  to  avoid  aliasing  effects  the  sampling  rate  must  be  chosen 
in  excess  of  2(10^)fQ,  say  f  =  (lO^Jf^.  A  million  samples  of  data 
for  both  input  and  output  are  thus  produced.  Apart  from  the  diffi¬ 
culties  of  storing  this  staggering  amount  of  data  and  the  impracti¬ 
cability  of  processing  them,  serious  numerical  difficulties  also 
arise  from  this  simple  minded  approach  to  identification;  for  in¬ 
stance,  the  low  frequency  poles  cannot  be  represented  in  z-domain 

accurately  even  with  a  64-bit  computer  word.  A  possible  remedy  is 

3 

to  break  the  problem  into  two  or  three  smallband  problems.  The 
network  dynamics  can  be  identified  for  each  of  these,  and  this  in¬ 
formation  can  be  used  to  estimate  the  wideband  transfer  function. 

A  mult  Irate  sampling  approach  to  identification  of  wideband 
transmittances  is  presented  in  this  section.  It  permits  efficient 

3 

A  frequency  band  of  less  than  two  decades  will  be  termed  as  small¬ 
band. 


34 


determination  of  the  transfer  function  of  a  four-to-five  decade 
bandwidth  system  from  transient  tests.  Clever  selection  of  sam¬ 
pling  rates  and  exciting  inputs  reduces  the  wideband  problem  into 
three,  simpler  smallband  problems.  Each  smallband  problem  encom¬ 
passes  only  one-to-two  decades  of  bandwidth.  The  three  transfer 
functions  H^(s),  H^(s) ,  and  H^(s)  are  easily  identified  via  the  pen- 
eil-of-functions  method,  and  then  adjoined  to  build  the  wideband 
transfer  function  estimate  H(s).  The  technique  is  demonstrated  by 
simple  illustrative  examples  and  a  realistic  UF  amplifier  example. 
Frequency  regions  (sub-bands) 

The  concept  of  small  band  descriptions  begins  by  splitting  the 
wideband  region  into  three  regions.  As  shown  in  Fig.  10,  these  re¬ 
gions  may  be  termed  as  low-frequency  band ,  medium  frequency  band, 

4 

ana  the  high-frequency  band;  in  short.  LF,  >1F,  and  HF.  These  re¬ 
gions  may  be  chosen  covering  approximately  equal  ranges  on  the  loga¬ 
rithmic  scale.  Denote  the*  band  edges  as  1  .  f  ,  f  ,  and  •  ,y  and  the 
respective  mid- region  frequenei.  s  as  • ,  ,  f  and  t.  .  The  latter  mar 
be  —  although  not  necessariJv,  cue  ,eti  as  the  geometric  means  of  the 
sand  edge  frequencies. 

By  design,  Lite  following  inequalities  hold 


o.if,  • 

f  •- 
o 

fL  '  fl 

10fL 

(40a) 

fl 

’n  '  '  > 

10fM 

(40b) 

1  t 

‘  ll 

t2  ' 

Si  ■  fj  ' 

l0tii 

( lUc) 

s 

in  some  cases  where  prior  knowledge  of  the  approximate  frequency 
characteristic  of  the  network  is  available,  it  may  be  more  appro¬ 
priate  to  choose,  the  regions  as  LF,  LMF  and  HF,  or  as  LF,  MHF  and 
HF. 
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Sampling  Rates 


For  Che  three  small  band  problems  the  sampling  intervals  are 
5 


chosen  as 


Al  -  1/lOOq 


A„  -  l/100f 

M  M 


Ah  =  l/100fH 


(Ala) 

(Alb) 

(Alc) 


The  sampling  rates  are  of  course  the  reciprocals  of  these  numbers. 

By  this  choice  —  and  in  view  of  (AO),  the  sampling  rates  become  at 
least  ten  times  the  highest  frequency  of  interest  in  the  respective 
bunds.  If  the  system  input  are  selected  so  as  to  excite  frequencies 
only  in  one  of  the  bands,  then  by  using  the  prescribed  sampling  rate 
aliasing  effects  wouid  be  avoided. 

Now,  if  K,  the  number  of  samples  used  for  the  identification 
procedure,  is  taken  as  1000,  tin  length  of  the  record  would  be  2 it 


times  tlie  longest  time  constant  of  tile  band  under  consideration  (for 


10 


example,  with  this  choice,  of  K  for  the  LF  case,  Tj  =  1000A^  =  —  = 
l 

—  2irr  ) .  Such  record  lengths  are  considered  adequate  for  practical 
o 


identification  of  low  edge  corner  frequencies,  and  storing  and  pro¬ 
cessing  1000  samples  of  data  is  well  within  today's  minicomputer 
capability . 

Inputs 

The  key  to  the  conversion  of  the  wideband  problem  to  three  small 
band  problems  is  the  careful  selection  of  inputs  which  excite 


These  are  conservative  values  in  anticipation  of  500  to  2000  data 
points.  Larger  values  up  to  five  times,  and  accordingly  fewer  data 
points,  may  he  used  with  some  caution. 
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frequencies  essentially  limited  to  one  of  the  bands.  At  first  this 
would  seem  to  pose  no  real  difficulty,  for  we  can  choose  a  narrow- 
band  signal  for  the  test.  However,  a  little  thought  would  reveal 
that  testing  with  very  narrowband  signals  would  be  in  direct  con¬ 
flict  with  the  basic  philosophy  of  system  identification,  which  is 
broadband  modeling  with  transient  tests.  Therefore  a  judicious  com¬ 
promise  must  be  made  between  these  conflicting  requirements. 

The  following  inputs  are  suggested  as  a  rough  guide.  Experi¬ 
mentation  and  experience  leads  to  a  much  richer  variety  of  signals 
which  meet  the  above  compromise  strategy.  Two  different  considera¬ 
tions  have  been  kept  in  mind  in  the  selection  of  these  inputs:  the 
spectral  requirement  stated  above  and,  equally  important,  easy 
realizability  in  the  laboratory, 
a)  LF  Input  - 

For  the  low  frequency  band  the  input  selected  is  a  triangular 
pulse,  either  a  full  cycle  TR  _(t)  or^  a  half  cycle  TR  (t)  (see 
Fig.  11).  In  either  case,  the  total  duration  of  the  pulse  is  taken 
to  be  T^/2  and  the  pulse  is  followed  by  zero  input  for  the  remainder 
of  the  time,  i.e.,  from  T  /2  to  T  .  The  magnitude  spectra  of  these 
inputs  can  be  shown  to  be 


SimrfTL/8 
nfT  /  8 

L 


2 

|Sin2irfTL/8| 


(42) 


I TR+ ( f )  I 


T  SirnrfT  /4 

_ L  _ la 

4A  irf Tt  /4 

L  la 


(43) 


where,  keeping  (15a)  in  mind,  T  =  K/lOOf  . 

la  la 


A™*-"1 

T.  0 

TR+(t> 

y  time 

L  u 

time 

(a)  (b) 

Figure  ll.  input  waveforms  for  LF  tests. 


For  networks  which  pass  d.c.,  TR,  (t),  i.e.,  a  full  cycle  triangular 

'  9~ 

pulse,  is  recommended;  this  reduces  the  predominance  of  a  d.c  component 
in  the  network  response. 
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Unit  peak  values  for  the  pulses  have  been  assumed.  The  amplitude  spectra 
are  tabulated^in  Tables  1  and  2 

TABLE  1 

Magnitude  Spectra  of  +,-  triangular  pulse  (  ) 


f/fL 

K=20 

|tr  (f)| 

K=200 

K=1000 

K=2000 

0.01 

-56.1  dB 

-36.1  dB 

-22.1  dB 

-16.1  dB 

0. 1 

-36.1 

-16.1 

-3.5 

-1.8 

0.5 

-22.1 

-3.5 

-16.1 

-29.8 

1.0 

-16. 1 

-1.8 

-29.8 

— OO 

2.0 

-10.3 

—CO 

—CO 

— oo 

10.0 

-1.8 

—CO 

—CO 

—00 

max  in 

band 

-1.8 

-1.8 

-1.8 

-1.8 

Note : 

zero  dB  corresponds 

to  a  magnitude  of  2b/A  =K/4 

Li 

TABLE  2 

Magnitude 

spectra  of 

a  +  triangular  pulse  (  /!_  ) 

1 tr+ ( f ) | 

f/fL 

K=20 

K=200 

K=1000 

K=2000 

0.01 

-0.0  dB 

-0.0  dB 

-0.0  dB 

-0.1  dB 

0. 1 

-0.0 

-0.1 

-1.8 

-7.8 

0.5 

-0.0 

-1.8 

-29.8 

-35.8 

1.0 

-0. 1 

-7.8 

-35.8 

— OO 

2.0 

-0.3 

—00 

— OO 

— OO 

10.0 

-7.8 

— oo 

— OO 

— oo 

max  in 

band 

0.0 

-0. 1 

-1.8 

-7.8 

Note:  zero  dB  corresponds  to  a  magnitude  of  2b/A  =K/4 

Lt 


Minus  infinity  is  used  whenever  iho  spectral  amplitude  is  less  than 
-200  dB  below  reference  level. 
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It  is  clear  from  Tables  1  and  2  that  the  spectra  of  these  pulse 
inputs  diminishes  to  -30  dB  or  more  (below  in-band  maxima)  at  the 

g 

LF-MF  boundary,  provided  N  is  chosen  greater  than  or  equal  to  200 
This  insures  that  the  frequencies  in  the  MF  region  are  not  excited  by 
application  of  these  inputs.  A  possible  exception  is  the  case  where 
there  is  a  sharp  resonant  peak  in  the  MF  band,  particularly  at  the 
LF-MF  boundary.  However,  the  presence  of  such  a  peak  is  generally 
known  before  hand;  such  a  resonant  component  in  the  output  can  be 
filtered  before  performing  identification  on  the  LF  test  data. 

MF  Input  - 

For  the  medium  frequency  band  the  input  selected  is  an  oscillatory 
pulse,  modulated  either  by  a  decaying  exponential  0EX(t)  or  by  a  dimini¬ 
shing  one-quarter-cycle  triangular  waveOT;>(t).  In  either  case,  the 
frequency  of  oscillation  is  taken  to  be  f^,  the  center  frequency  of 
the  band.  The  duration  of  the  oscillation  is  taken  to  be  T^/ 2  (see  Fig.  12) , 
followed  by  zero  input  for  the  remainder  of  the  time,  i.e.,  from  T../2 
to  T^.  In  presenting  the  spectral  analysis  below  it  is  assumed  that 
the  on-set  of  the  pulse  begins  with  the  maxima  of  the  oscillation,  i.e., 

the  pulse  is  triggered  at  its  maximum  value.  Thus  u(t)=  m(t)  Cos  2iTf  t 

M 

where  m(t)  denotes  the  modulating  envelope.  The  spectra  of  these  inputs 
can  be  shown  to  be  as  follows: 


|  OEX ( f ) |  =  -y-  |M(f+fM)  +  M(f-fM)|  , 


1  r  -Ca+ju»)T/2 

M(f)  =  ri  - 6  ] 


Kf)  |  =  4“  lM(f+fM)  +  M(f_fM)l  • 


M(f)  =  ~r~~  [  1 
jw 


SinrrT^f/2 

ttT  f/ 2  (CoslTTMf/2  -  JSinirT  f/2)] 
M 


sIf  the  sampling  interval  were  chosen  five  times  the  value  suggested 
in  (Ala)  the  magnitude  spectrum  diminishes  to  -30  dB  at  the  LF-MF 
boundary  even  for  N=20. 

AO 


ST 


where  u)=2Trf,  'a'  is  the  inverse  time-constant  associated  with  the  expo¬ 
nential  decay  and,  keeping  (41b)  in  mind,  T  =K/100fw.  Unit  peak  values 

MM 

have  been  assumed. 


Figure!!.  Input  waveforms  for  MF  test. 

In  order  to  delineate  the  spectral  characteristics  of  the  input 

OEX(t) ,  three  different  values  of  'a'  will  be  considered:  a=0,  a=2/Tw, 

M 

and  a=4/T^.  The  corresponding  waveshapes  will  be  denoted  as  OEXp(t), 


OEX1 (t)  an 
OEXx(t),  0 

f/fM 

id  0EX2(t),  respectively.  The  amplitude  spectra  of  OEX^ft), 
iEX2(t)  and  OTR(t)  are  tabulated  in  Tables  3  to  6  respectively 

TABLE  3 

Magnitude  Spectra  of  an  Oscillatory  pulse 

K  =  20 

|0EXq 

K  =  200 

(Ol 

K  =  1000 

K  =  2000 

0.010 

0.1  dB 

-74.0  dB 

-74.0  dB 

-74.1  dB 

0.100 

0.1 

-34.0 

1 

OJ 

•^J 

00 

—  OO 

0.500 

0.1 

-7.4 

-21.4 

—  00 

0.909 

0.0 

-0.5 

-3.6 

-20.5 

0.990 

0.0 

-0.0 

-0.1 

-0.2 

1.000 

0.0 

0.0 

0.0 

0.0 

1.010 

-0.0 

0.0 

-0.0 

-0.1 

1.100 

-0.0 

0.3 

-3.5 

—  OO 

2.000 

-0.4 

— oo 

—  OO 

—  00 

10.000 

-29.1 

—  03 

—  oo 

—  OO 

Note:  zero  dB  corresponds  to  the  resonant  peak  at  f  .  Minus 

M 

infinity  is  used  whenever  the  spectral  amplitude  is  less 
than  -200dB  below  reference  level. 

41 


•  > 


i 


TABLE  4 


Magnitude  spectra  of  an  Oscillatory  exp.  pulse 


f/fM 

K  =  20 

|OEX1(f) 

K  =  200 

. 

i 

K  =  1000 

K  =  2000 

0.010 

-1.7  dB 

-26.2  dB 

-53.1  dB 

-63.1  dB 

0.100 

-1.6 

-23.4 

-36.7 

-49.7 

0.500 

-1.1 

-6.9 

-20.7 

-33.5 

0.909 

-0.2 

-0.5 

-3.4 

-14.6 

0.990 

-0.0 

-0.0 

-0.1 

-0. 2 

1.000 

0.0 

0.0 

0.0 

0.0 

1.010 

0.0 

0.0 

0.0 

-0.1 

1.100 

0.2 

0.3 

-3.2 

-15.7 

2.000 

1.2 

-13.6 

-27.4 

-33.5 

10.000 

-12.9 

-29.9 

-43.8 

-49.9 

Note:  zeru  dB  corresponds  to  the  resonant  peak  at  f 

TABLE  5 

Magnitude  spectra  of  an  Oscillatory  exp.  pulse 
i  OEX.;  (  f )  | 


l/fM 

K  =  20 

K  =  200 

K  =  1000 

K  =  2000 

* 

0.0  L0 

-0.0  dB 

-15.0  dB 

-41.7  dB 

-53.2  dB 

0.100 

-0.0 

- 14 . 2 

1 

C<0 

■O 

o 

-43.4 

O 

V 

0.500 

-0.0 

-5.6 

-19.1 

-27.4 

> 

0.909 

-0.0 

-0.4 

-2.9 

-9.8 

0.990 

-0.0 

-0.0 

-0.1 

-0.2 

t 

1 . 000 

0.0 

0.0 

0.0 

o 

o 

*1 

1 . 0 1  0 

0.0 

0.0 

0.0 

1 

c 

f— 1 

\ 

*4 

1.  100 

0.0 

0.2 

-2.6 

1 

O 

o 

2.000 

-0.1 

-.8.  1 

-21.4 

1 

ro 

-IS 

< 

1 0 . 000 

-9.3 

-24.1 

-37.8 

1 

CO 

00 

» 


Note:  zero  dH  corresponds  to  the  resonant  peak  at 

42 


TABLE  6 


Magnitude  Spectra  of  an  Oscillatory  triangular  pulse 


f/fM 

K  =  20 

| 0TR(f ) | 

K  =  200 

K  =  1000 

K  =  2000 

0.010 

0.1  dB 

-38.0  dB 

-51.9  dB 

-58.1  dB 

0.100 

0.1 

-18.0 

-36.3 

-43.8 

0.500 

0.1 

-4.3 

-21.2 

-27.4 

0.909 

0.0 

-0.3 

-2.3 

-8.7 

0.990 

0.0 

-0.0 

-0.0 

-0.1 

1.000 

0.0 

0.0 

0.0 

0.0 

1.010 

-0.0 

0.0 

0.0 

-0  1 

1.100 

-0.0 

0.2 

-2.1 

-9.5 

2.000 

-0.3 

-7.6 

-21.4 

-27.4 

10.000 

-9.3 

-23.9 

-37.8 

-43.8 

Note:  zero  dB  corresponds  to  the  resonant  peak  at 

It  is  clear  from  Tables  3  and  4  that  the  spectra  of  pulses  OEX^t)  and  OEX^ 

(t)  diminish  to  -35  dB  or  more  at  the  Mf-LF  and  MF-HF  boundaries  provided 

q 

N  is  chosen  greater  than  or  equal  to  200  .  This  insures  that  the  frequencies 

in  the  LF  and  HF  regions  are  not  excited  significantly  by  application  of 
these  inputs.  Tables  5  and  6  show  that  0EX2  and  OTR  spectra  diminish  only 
to  -15  dB  at  these  boundaries;  these  pulses  are  useful  in  the  initial 
stages  of  testing,  or  when  the  network's  corner  frequencies  are  spread 
over  the  band. 

HF  Input  - 

The  inputs  used  for  the  MF  band  are  equally  useful  for  the  HF  band  with 

T  replaced  by  T„.  Tables  3  to  6  also  hold  with  f  replaced  by  f  . 

M  H  M  rt 


9  If  the  sampling  interval  is  chosen  five  times  the  value  suggested  in  (41b), 
the  magnitude  spectrum  diminishes  to  -30  dB  at  the  boundaries  even  for  N=20. 
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Small  band  Identification 

Input-output  data  obtained  from  smallband  tests  can  be  analyzed  by 
use  of  the  Fortran  program  IGRAM  [5]  and  the  s-domain  smallband  transfer 
functions  obtained  therefrom.  The  program,  however,  requires  that  the 
transfer  function  order  (degree  of  the  denominator  polynomial)  be  speci¬ 
fied.  If  the  order  of  a  smallband  transfer  function  is  known  from  cir¬ 
cuit  considerations,  then  the  identification  is  performed  for  this  order 
and  for  at  least  one  order  higher  and  lower.  For  example,  if  the  LF  band 
behavior  is  expected  to  be  of  order  4,  then  identification  should  be  per¬ 
formed  for  n=3,  4  and  5.  The  lowest  order  model  yielding  satisfactory 
fractional  error  (see  page  15)  should  be  accepted  as  the  model  for  that 
smallband.  If,  on  the  other  hand,  the  smallband  order  is  not  known,  then 
an  upward  modeling  strategy  must  be  adopted.  Starting  from  an  initial 
order,  a  low  guess,  increasingly  higher  orders  are  attempted  until  the 
fractional  error  in  identification  turns  out  to  be  acceptably  small. 

Thus,  the  smallband  transfer  functions  H^(s),  H^(s)  and  H^(s)  be¬ 
come  available.  From  these  the  overlapping  critical  frequencies,  or 
ideally  speaking  common  critical  frequencies,  are  carefully  isolated. 

This  isolation  of  common  critical  frequencies  is  useful  in  the  next,  and 
final,  step  in  wideband  identification. 

Adjoined  Wideband  Transfer  Function 

The  transfer  functions  obtained  from  the  smallband  tests  must  he  ad¬ 
joined  to  form  the  overall  wideband  transfer  function.  For  convenience  we 
will  drop  the  hat  (carat)  on  the  identified  TFs,  smallband  or  wideband. 

The  reader  must,  however,  bear  this  in  mind. 

In  the  notation  to  follow  we  will  use  C  iO  denote  gain  constant;  and 
H  with  suitable  subscripts  to  denote  transfer  functions,  which  are  assumed 
to  be  in  the  Bode  canonical  form 

s'"(l  +  s/z.)  ....  (1  +  s/z  ) 

C _ — _ I - SL 

(1  +  s/p^  ....  (1  +  s/pn) 

Ttie  first  subscript  (on  H)  refers  to  the  test  from  which  the  transfer 
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''' 


(46) 

(v  =  0,  positive 
or  negative  integer  ) 


•  s 


T 


function  is  obtained;  the  second  subscript,  if  any,  denotes  the  band  with 
which  the  critical  frequencies  (poles  and/or  zeroes)  are  shared.  For 
example,  H^(s)  denotes  the  part  of  H^(s)  whose  critical  frequencies  are 
shared  with  H  (s) .  Ideally  speaking,  of  course,  H^ffl(s)  = 

Thus  we  have 


Vs)  ■  CL  V<a)  Vi<s)  <47-» 

VS>  ■  CH  W0  VS)  "»)<=>  <47b> 

vs>  - c;  «,«<=>  v<s)  <47c> 


Critical  Frequency  Adjustment  - 


In  practice  the  overlapping  critical  frequencies  in  two  small- 
band  tests  will  not  turn  out  to  be  identical  through  the  corresponding 
identifications.  For  example, a  critical  frequency  s|  (pole,  or  zero) 
common  to  H^(s)  and  H^s),  may  be  identified  as  s^_  in  the  MF  identifica¬ 
tion  and  as  s^+  in  the  HF  identification.  We  will  adjust  both  of  them  to 
a  common  value  given  by  their  geometric  mean 


s.  =  J 


Si+ 


(48) 


Assume  that  this  process  has  been  performed  on  and  H^,  and  like¬ 
wise  on  and  H^.  In  order  to  avoid  unduly  complicated  notation  we 
will  let  the  original  symbols  denote  these  adjusted  transfer  functions, 
so  that  now 


V>  ■  V,(s> 


(49a) 

(49b) 


Other  adjustments  include  setting  s^  to  0  when  it  turns  out  to  be  well 
below  2nf  ,  where  f  denotes  the  left  boundary  of  the  frequency  band,  but 
is  known  to  be  zero  from  circuit  considerations;  care  should  be  taken  in 
this  case  to  let  s^(l  +  s/s^)  to  s.  Thus,  in  the  canonical  form  of  (46) 
the  gain  should  be  divided  by  s.  when  the  term  (1  +  s/s^)  is  replaced  by 
s.  Another  case  of  adjustment  occurs  when  s^  turns  out  to  be  much  larger 
than  2irf+,  where  f+  denotes  the  right  boundary  of  the  frequency  band;  it 
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is  Chen  useful  to  set  the  term  (1  +  s/s^)  to  just  1.  These  two  types  of 
adjustments  occur  in  Example  8  on  page  47. 

Gain  Adjustment 


To  obtain  equalization  at  the  boundaries  of  the  frequency  bands  the 


gains  are  adjusted  as  follows. 

„  IVVWI1'2 


°L  -  CI 


vvvv 


(50a) 


r  r,  vwv 


(50b) 


„  ww; 

“  VVW 


where 

=  j2irf 

A2  =  j2TrE2 

Recall  that  f  is  the  boundary  between  the  LF  and  MF  regions,  and  f^  is 
the  boundary  between  the  MF  and  HF  regions  (see  Fig.  10,  page  36). 

As  stated  earlier,  the  purpose  of  this  gain  adjustment  is  to  minimize 
gain  discontinuity  at  the  boundaries.  However,  phase  mismatch  may  still 
exist  at  these  boundaries  for  the  redefined  smallband  transfer  functions. 
These  transfer  functions  are 


Ht  =  c,  H, ,  (s)H  (s) 
L  L  LL  LM 


(51a) 


VS)  ‘  CM  HML(3)lWS)lWS) 


V“>  ■  c»  ^HM^S^HH^S ^ 


(51c) 


Wideband  Transfer  Function  - 

The  wideband  transfer  function  is  taken  to  be 

•  c  \l<s) Vs)Ws)hhm(s)hhh(s)  <52) 

For  certain  wideband  networks  only  two  smallbands,  LF  and  HF  (with  boun¬ 
dary  frequencies  f  ,  f^  and  ^2)9  might  be  necessary.  In  such  cases  only 
(50a)  and  (50c)  are  needed  with  H^  deleted  and  X ^ = A 2 = j 2nf ^  ;  likewise,  onl 
(51a)  and  (51c)  are  needed  wherein  the  subscript  M  is  replaced  by  the  sub 
script  H. 
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where  the  constant  C  is  selected  to  match  the  gain  of  one  of  the  smallband 
transfer  functions  at  a  chosen  frequency,  perhaps  H^s)  at  the  midband 
frequency  f^. 

Example  8  - 


Consider  an  a.c.  coupled  network  believed  to  have  frequencies  of  in¬ 
terest  from  0.02  MHz  to  50  MHz,  thus  encompassing  3.3  decades.  The  para¬ 
meters  of  the  system  are  given  in  Table  7.  In  this  case  it  is  adequate  to 
break  up  the  frequency  region  into  two  smallbands  as  follows: 

LF  f  =  0.02  MHz  to  f.  =1.0  MHz  fT  =0.1  MHz  (53a) 

o  1  T 


HF 


f1  =  1.0  MHz  to  f2  =  50.0  MHz 


f  =  0.1  MHz 

L 

f  =  10.0  MHz 
H 


(53b) 


The  inequalities  in  (40)  are  clearly  satisfied.  From  (41)  and  Tables  1 
and  3  it  appears  reasonable  to  choose^ 

=  200 


A^  =  0.1  us 


Tl  =  20  ps 


A  =  0.001  ps 
n 


"h  "  200 


T„  =  0.2  ps 


and  the  inputs  u  (t)  =  TR  (t),  u  (t)  =  0EX  (t).  We,  however,  select 
L  •  > —  H  U 


A^  =  0.05  ps 


*L 


=  200 


Ah  =  0.002  ps  Kh  =  100 


T  =  10  ps 
L/ 

Th  =  0.2  ps 


(54a) 


and  define  the  test  pulses  explicitly  as  follows 

u^(t)  =  f  One  complete  cycle  triangular  wave  over  0  to  5  sec.  (55a) 
[0  level  over  5  to  10  sec. 


UH(t)  '  Cos(-— 2gy-  t),  0  <  t  <  0.1  (55b) 

0  0.1  <  t  <  0.2 

To  simulate  the  LF  test  response,  the  system  function  ^true(z)  corres¬ 
ponding  to  A'  =  0.005  psec.  is  excited  by  u  (k)  of  (55a);  the  response  is 

L  L 

then  resampled  at  l/10th  rate  (i.e.,  every  10th  output  sample  is  picked 
up).  In  a  laboratory  test  this  artifice  of  using  a  high  sampling  rate 
H(z)  to  preserve  the  integrity  of  the  network  response,  and  then  resampling 
the  output,  is  of  course  not  necessary.  The  network  output  can  be  sampled 
directly  at  the  desired  rate  1/A  . 

^  Recall,  A  denotes  sampling  interval,  K  the  total  number  of  samples  and 
T  the  total  duration  of  the  test.  Of  course  T  =  KA. 
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Table  7 

A  Wideband  System:  Example  8 


H  (S)  =  _ _ 

true  }  (s+106)(s+108) 


A^  =  0.00 5ps 


H  (z)  = 

true,  i 


0.003924z~1q-z~1) 
1-1 . 601543z-1+0. 603505z 


A  =  0.002ps 
n 

H  (z)  =  0.0018108z~1Cl-z~1) 

tru6,AH  l-1.8167327z_1+0.81709493z' 


The  results  of  identification  from  IGRAM  using  the  method  of  Subsec¬ 
tion  2.1  are  given  below. 

LF  Test  - 


Both  first  and  second  order  identifications  were  performed.  Since 
the  first  order  model  (predicted  response)  has  a  fit  error  of  v  *  0.8E-4, 
it  is  decided  the  LF  behavior  is  first  order.  For  n=l  the  Gram  matrix, 
the  square-roots  of  the  diagonal  cofactors,  the  z-domain  model  and  the 
corresponding  s-domain  model  (using  pulse-invariant  conversion  [5],  pages 
80-82)  are  given  below 

Hj. 00074285  0.00037153  -0.075250  5.1001 

1.03620000  -5.198800  -354.5100 

33.336000  1266.7000 

|  346090.00 


det  F  =  0.128 


The  values  of  /D^/D^  are 


0.050249 


0.010103  -0. 7329E-6 


yz) 


=  0.0096187 (1-z 
l-0.95215z“] 


Vs) 


0.0101 (s-0 . 14 (10 
(s+0.981(106)) 


1 . 0297 (10~8) s 
(s/0.981(106)+l) 
48 


(Frequency 

adjusted) 


p*  - 
*  & 


HF  test 


As  in  the  LF  case,  here  also  a  first  order  model  is  found  adequate 
producing  a  fit  error  (fractional  energy  error)  v  =  0.40E-4.  For  n=l  the 
Gram  matrix,  the  square-roots  of  the  diagonal  cofactors,  the  z-domain 
model  and  the  corresponding  z-domain  model  (using  pulse-invariant  conver¬ 
sion)  are  given  below. 


0.00090304  0.00045161 

0.01842400 


det  F  =  0.22E-4 


The  values  of  /D^ /D^  are 


-0.073612 

0.449400 

25.000000 


-0.52301 

-1.61760 

12.50000 

410.38000 


0.22817 


-0.0023153 


0.0022408 


Hh(z) 


-0.00006  +  0-00188517Z 
1  -  0.814217z-1 


Hh(s) 


0.0023s  +  1.009(10  )  _ 
s  +  102 . 76 (106) 


0.00982 

s/1.028(108)  +  1 


(Frequency 

adjusted) 


Gain  Adjustment  - 

At  s  =  j 2n f ^ ,  where  f^lMHz,  the  gains  of  the  LF  and  HF  transfer 

functions  turn  out  to  be  0.00998  and  0.00981  respectively.  The  adjusted 

gain  constants  (using  (50)  and  (51))  are  C=l. 021(10  8)  and  C  =0. 00990. 

L  H. 

Wideband  Transfer  Function  - 

The  wideband  transfer  function  is 


H (s)  =  C 


(s/0.981 (106)  +  1) (s/1.028(108)  +  1) 


where  C=l. 023(10  )  is  obtained  from  gain  matching  at  s=j2nf^;  f^=lMHz. 

Comparison  - 

The  Bode  plots  of  H(s)  and  Htrue(s)  (of  Table  7)  are  compared  in  Fig. 

It  appears  that  satisfactory  wideband  identification  has  been  achieved. 

Remark 


The  procedure  of  adjoining  the  smallband  transfer  functions  can  of 


course  be  programmed. 


As  a  second  example  of  wideband  identification  consider  the  RF  ampli¬ 
fier  (Fig.  14)  of  reference  [5].  The  frequency  regions  are 


f  =0.002  MHz  to  f.=0.1  MHz 
o  1 

f1=0.1  MHz  to  f 2=10  MHz 
f 2=10  MHz  to  f3=1000  MHz 


f=0.01  MHz 
L 

f  =1.0  MHz 
M 

f  =100  MHz 
H 


The  smallband  transfer  functions,  identified  from  LF  MF  HF  tests  through 
I GRAM  [5]  are 


Vs)  = 


=  -20.125 


(s-0. 0015  (106) )  (s-t-0. 0012  (106)) 
(s+0.034(106))(s+0.075(106)) 


-0.7892(10  8) 


[s/0.034 (106)+1H s/0. 075 (106)+1] 


H  (s)  =  -^20(106)  - - - 7 — 

(s+24. 92(10°)) 


-20. 3  5 


(s/25 . 31 (10°)+1) 


=  2.79(10  ) 


_ (s-19060 (10°) ) _ 

(s+25. 7 (106) ) (s+1140(106) ) 


18.432 


_ (-s/19606 (10  )+l) _ 

[ s/25. 31 (106)+1] [s/1140 (106)+1] 


The  second  step  of  each  of  the  above  is  obtained  after  frequency  adjustment 
as  outlined  on  page  45. 

Gain  Adjustment  - 

At  f^=0.1  MHz  the  gains  of  the  LF  and  MF  transfer  functions  turn  out 

to  be  19.9536  and  20.5437  respectively.  At  f2=10  MHz  the  gains  of  the  MF 

and  HF  transfer  functions  are  computed  to  be  7.6776  and  6.8759  respectively. 

The  adjusted  gain  constants  (using  (50)  and  (51))  are  C  =(0.9602)C'=0. 7578, 

L<  L 

CM=(0. 9327)0^=19.17,  and  CH=(1. 0414)0^=19.19. 

Wideband  Transfer  Function  - 


The  transfer  function  of  the  network  is  estimated  as 


H(s)  =  -C 


_ s2[-s/ 19060 (106)+1] _ 

[s/0.034(106)+l] [s/0.075(106)+l] [s/25.31(106)+l][s/1140(106)+l] 


(a)  Schematic 


(b)  Equivalent  circuit 


C^O.OI  Ff 

C  =  5  pF  g^=4  mmho 

g^=A0  mmho 

Z  =  1  Kohm 
(real ) 

C2=  50  pF 

C^=0.01  g2=1  mmho 

g^=0.5  mmho 

8 (TO7)  s2(s-8000(106)) 

( s+0 .03 3  ( 1  f'6 ) )  ( s+0 .08  (106) )  ( s+25 . 2  ( 1 06 ) )  ( s+ 1205 .1(1 oS ) 


-  7-983 (10~9) [ s/8000 .0(10°)  +  1  ] 


I  s/0.033 ( loS  +  1 ] [ s/0. 08{ 106)  +  I][s/25.2(106)  +  1  ]"[s/1205. 1  ( 1 06)  +  1] 


Fig.  14  A  common-emitter  wideband  amplifier 


TT 


_Q 

where  C=8. 0596(10  )  is  obtained  from  gain  matching  with  H^s)  at 

f  =1  MHz. 

M 

Comparison  - 

The  Bode  plots  of  H(s)  and  ^true(s)  (of  Fig*  14)  are  compared  in  Fig. 
15.  It  appears  that  satisfactory  wideband  identification  has  been  achieved 
from  smallband  time  domain  tests. 
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•  s 


Comparison  of  the  magnitude  (Bode)  plots  of  the  identified 

transfer  function  and  the  network  function 
of  an  RF  amplifier 
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oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


STOZ 


GIVEN  the  CONTINUOUS  CESCFIPTICN,  PROGRAM  CO“PuTEi  Tnf 
EQUIVAlENT  OISCPETC  o OM A  IN  DESCRIPTION  0C  <•  LINEAR 
OYNAMIC  SYSTEM 

STOZ  GENERATES  H(Z>  ANC  TNF  CG  r  RESPOND  IN  G  Cl  ff  c°  -'IJCE 
EQUATION  FROM  THE  TRANSFER  FUNCTION  HCS) 

THE  INPUT  ARRAYS  A  ANC  ft  A  FE  FILLED  ACCORrlJC  TO  THE 
CIFFFRNTIAL  EQUATION 

ft  (1>»  Y(T)  ♦P(2»  *r>Cl,  Y(  T  »  I  ♦  ..  .  *e(N»l)  *0  <N,Y  «  T  I  > 

-A  (1 )  *U(T)-A  (2>*D(1,  U  (T)  >-.  .  AIN*  1 »  *  C  (N,U  (T  *) 
WHERE  0(M,F(T)>  =  THE  NTH  TINE  OERIV  AT I V  E  OF  FUNCTION,  F 

P (N*  1 )  HUST  EQUAL  1 

RETURNS  ARRAYS  A  ANC  e  CONTAINING  THE  EQUIVALE'T  QJSf.AF- 
CESCRIRTION  STORED  ACCORDING  TC  THE  OIFFER£nC£  rCjAIION 

0(l)*Y(K)4B(2)»Y(X-:M...*ft(N41)»Y(K-N> 

-A  (1)  *tl(K| -A  (2)»L(K-1) -•  ..-A  <N +1 )  *  L  (  K -N)  =  0 

8(1)  ALWAYS  EQUALS  1 

THE  ROIES  OF  T«c  CONTINUOUS  DOMAIN  HUST  BE  DISTINCT  A- C 
NON-ZERO  FOR  The  TRANSFORMATION  TO  BE  VALID 


DATA  CARO  SET  PREPARATION 
N  =  ORDER  CF  SYSTEM 

n  (MAXIMUM)  =  ONE  LESS  than  The  OIKENSION  SUSSC^I0* 
iht h o  i  o  for  Twr  impulse  invariant  description  . 

=  1  FOR  ThE  PULSE  INVARIANT  DESCRIPTION  . 

=  2  FOR  ThF  TRAFf7CIr'AL  INVARIANT  DE  S.C  SI*  7  TON  . 

=  3  FOR  THE  LOGRITHMIC  TRANSFORM  QEC  DR  1PT I  Ctl  . 

IPOLZ  =  1  IF  POL  rS  ANC  ZEROES  A ®E  READ 

MUST  BE  COMFLEX  (RE A L . IM AGIN AR Y ) 

NEGATIVES  OF  PCLES  AND  ZFRCS  READ!  IF, 

FOR  A  POL  Z  CF  (S  +  2),  INPUT  +  2.E  ♦  C  .  j 

(2F10.C  PEP  FOLZ,  A  FOLZS  PER  CARD) 

=  C  IF  DENOMINATOR  AND  NUMERATOR  ARt  pf AT 
IN  POLXNOMIAL  FORM. 

COEFFICIENTS  ORDERED  FROM  LOW  TC  HIGH  DEGREE, 
DENOMINATOR  INFORMATION  ALWAYS  PEAC  FIRST. 
HIGHEST  ORDER  DENOMINATOR  COEFF  MUST  ftc  1.0 

NOTE  t  WHEN  ROLE-ZERO  CATA  IS  ENTERED  A  GAIN  CARD 
MUST  FOLLOW  THE  LAST  POLF  G..RC. 

IF  THERE  ARE  NO  ZFPOS,  USE  PLANK  CAPOS  IN 
THE  NORMAL  ZEROS  POSITIONS. 

CELTA  =  SAMPLING  INTERVAL 
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c 

c 

C  START  EACH  DATA  CAPO  SET  WITH  A  OESCRIPTICN  C ;<0. 

C  CONTAINING  UP  TO  51  CHARACTERS  COLS  2-52 

C  FIRST  DATA  CARD  CONTAINS* 

C  N,  IM-NC.  IPOLZ.IN  JF5  FCRMAT,  PLUS  DE  LTA  In  RlC.f  rc-'"A 

C  SECOND  GROUP  OF  DATA  CARES  IS  N  POLES  CR  f.  +  l  Pc n OM'n AT  C-' 

C  COEFFICIENTS.  AFTER  LAST  POLE  CARO,  JSt  A  C-aTn  CARD. 

C  LAST  GROLP  OF  DATA  CARDS  IS  N  ZEROS  (OR  FLANKS). 

C  OR  N  ♦  1  NUMERATOR  COEFFICIENTS  (BLANK'S  FC=  ZERO  0  5  F  F  S  * 

C  THE  OATA  FORMAT  FOR  EACH  OF  THE  SYSTEM  PARAMETER  CA^Dj 

C  IS  8F1C.C- 

C 

C  AS  MANY  SETS  OF  DATA  CARCS  MAY  BE  FUN  AS  EES!0", 

C  - - - - 

C 

C  ST 0 Z  MAIN  PROGRAM 

C 

REAL  3(20). »(2D  >, RR(2C > , K I (2C ) .CELT A, TEMP  (2  0> 

COMPLEX  C-  (21  >,  CA  (20)  ,C9  (20,  CAA(2i>  ,CA1(Z:  I  , 

1  TEM(20),C3N1,CON2,CONT 
DIMENSION  T IT  LE  (70) 

C 

C  RE  AO  TITLE  AND  ficsi  DATA  CARD 

C 

IOC  REAC(S ,B2  , 1 TItLE 

IF  (EOF  (5)  .N£.  C>  GO  TO  5939 
wRITE(fc,9lC) 

91  C  F  CRMAT  (6  (  /)  ) 

WRITE (6,920 )T ITLE 
REAC(5,92')TITLE 
92C  FCRMAT(ZCU) 

WRITt(6,973> 

93 C  FCRMAT(/,lX,71(*t*l( 

REAC(5,9<,3)  N,  IMTho.IRCLZ.  CELT  A 
9A0  FORMAT (313, Fit. C) 

WRITE (6,9^0 (N.IMTHD.IPOLZ.DELTA 
950  FORMAT  (//3X  ,*N  =•  ,  I<* .  5  X,  » IM.TmC  =  *  »Ih*5X  ,  * I»  3L  Z  =»,!-,  Rx, 

1  *  OELT  A  =  «,G17.U,//) 

NPl  =  NH 
NP  Z  =  N+  2 
NPNF1=NkN*1 
HPHPZ  =  H*ti*Z 

IF  ( IPOLZ.NE.l  )  GO  TO  300 
C 

C  RE  AO  PCLEE  ANO  ZSRCS 
C 

REAO(5,960)  (C  R  C I  >  ,  1=  1 ,  N> 

960  FORMAT  (  5F  lfl  .0  ) 

CALL  PClCCN(OR,  TEM,0,M 
00109  1=1, NP1 
109  B(II=TEm(I) 

C 

C  READ  GAIN  CARO 

R£AG<5,96  0>  P K 
C 

C  READ  ZEROS 

RE  AO  (5 , 96  0  )  (CA  (I)  ,1  =  1, N) 

CALL  POLCON (CA, TE M,0,N) 

DO  209  1=1, NP1 
209  A(I)=T£WU)*RK 

GO  TO  310 
C 

C  READ  DENOMINATOR  ANO  NUMERATOR  COEFFICIENTS 
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c 

30  C  RE  A  C  (5  »  96C  I  (  B  ( II  ,  1=  1 ,  NP  1 ) 

REAG(5»96Q)  <  A (  I)  ,  1=  1 ,  NP  1  > 

31 C  CONTINUE 

c 

C  PRINT  OrNOMINAT  0°  AND  NUMFPATOP  COEFFICIENTS 

C 

WRITF(6,9Z0> 

97 C  FCRM«T(*  S-OOMAIN  DENOMINATOR*) 

CALL  PRV£C(E,NP1I 
WRITE(6,91QI 

980  format  (*  s-oomain  numerator*) 

CALL  PRVECU.NPII 
C 

C  DETERMINE  ordfr  GF  NUMERATOR 

c 

NN  =  N 

DO  309  1=1,  NP1 
II=NP1*1-I 

IF  IKIII.NE.O  .0  1  GO  TO  40C 
309  NN=NN-1 
40  C  CONTINUE 

HRITE(6,993)  NN 

99  C  FORMAT!*  ORDER  CF  S-OOMAIN  NUMERATOR  =».I5.//> 

IF (NN.LT , 31  GO  TO  5079 
NNP1=NN*1 

c 

C  FACTOR  9EN0MINAT0R  TO  FIND  POLES 
C 

IF (IPOIZ.NE.O l  GO  TO  501 
CALI  PCLR'(e, TEMP,N,RR,RI,  IER) 

OO  4C9  1=1, N 

409  CR( II=CMPLX (OR(I) ,51(1)1 
50  C  WRITf(6,991> 

991  FORMAT (*  POLES  OF  T«F  S-CCMAIN*) 

CALL  P5CVFC (CP, M 

c 

IF(IMTMO.NE.T)  GO  TO  150  C 
C 

C  LOGRITHMIC  transform 

c 

100C  CONTINUE 

WRITE  (  6 , 9 0  G  • 

91 OC  F  CRMAT ( /,  *  LOGRITHMIC  TRANSFORM*) 

C 

C  WORK  ON  NUMERATOR 

c 

IF(NN.EO.G)  GO  TO  103C 

IF  (  IPOIZ.NE.O  )  GO  TO  1010 

CALL  PCLRT( A, TEMP,NN,RR,RI, TCP) 

00  1009  1=1, NN 

1009  CA<II  =  CMPLX(RR<I),RUI>I 

1010  HRITE(6,9?0C) 

9200  FORMAT ( »  ZEROS  IN  S  OONAIN*) 

CALL  P9CVEC(CA,NN) 

00  1029  1=1, NN 

1029  CA(II  =  CEXR(CA  (II*CELTA) 

IF(nn.EQ.N)  GO  TO  HOG 

1030  00  1039  I=NNP1,NP1 
CAA(II=U.GOO 

1039  CA(II=0.0 

1100  CONTINUE 
C 

C  NOW  THE  FIRST  NN  ENTRIES  OF  CA  CONTAIN  Th  r 


C  Z-OCMAIN  ZEROS  OF  THE  TSANSFFF  FUNCTION,  WHILE  7HE 

C  REMAINING  ENTRIES  AAE  ZEFCEC  CUT. 

c 

C  WORK  ON  DENOMINATOR 

C 

00  1129  1=1. N 

1129  CS<I>=CEXF«CRtI>*DELTA> 

C 

C  NON  C R  CONTAINS  THE  N  Z-COMAIN  POLES 

c 

C  FORM  NUMERATOR  AND  DENOMINATOR 
C  Z-OCMAIN  POLYNOMIALS 

c 

IF(NN.EO.C»  CAA(1>=1.C 

if<nn.ne.ci  call  pcstz(ca,caa.:.nn> 

CALL  PCSTZ <  CR  ,  C9,  0  »N  > 

C 

C  NON  CB  CONTAINS  THE  NN1  Z-DOMAIN  CENO“I NATO R  COEFFICIENTS. 

C  ANO  CAA  CONTAINS  THE  NN*  1  NUMERATOR  CC£FF IC IE  NT S . 

C 

C  ADJUST  DC  GAIN  CONSTANT 
C 

A1=A<1)/9<1> 

A2=1.3 

DO  1209  1=1, NN 
1209  A 2  =  A2* C  AA  ( I  +1  > 

92=1.0 

00  1219  1=1. N 
1219  92=a2*CP<T*l» 

F  A  C  =  A1 *  92/A2 
DO  1229  1  =  1  .NfJPl 
1229  CAA<I»=CAUI>»FAC 
C 

C  NOW  CAA  CONTAINS  THE  ADJUSTEC  Z-OOMAIM  N  Um  “R  kT  C°  COEFFICI-HT 

C  AND  FAC  CONTAINS  THE  GAIN  FACTOR  USED  Ff.R  th£  AOJ JSTf  ENT 

C 

GO  TO  500; 

c 

1500  CONTINUE 
C 

C  NON-LOC-RITHMI  C  transformations 
C 

C  ADJUST  F0=>  DIRECT  TRANSMISSION 
C  THIS  ROUTINE  REQUIRES  THAT  B(N»1I  =  1.0 
C 

CONT=  <0.5,0  .0 > 

IF<NN.LT.  it  GO  TO  151C 
CONT=A  (NR1» 

001509  1=1, N 

15  C  9  A«I>  =  AII>-CONT*BU) 

C 

C  FIND  NUMERATOR  CONSTANTS  FOR  PARTIAL  FRACTION  EXPANSION 
C 

1510  001529  1=1, N 

CON1=1.0 
CON  2=0 . 0 
001519  J=1,N 
C0N2=C0N2*Cfif I>NA(N-J+1» 

IFII-J)1512, 1519, 1512 
1512  CON1=CON1» (CR f I » -CR  t J) * 

1519  CONTINUE 
1529  CA ( I»=CON2/CON1 
WRITE<6, 93001 

93 OC  FORMAT <*  NUMERATOR  CONSTANTS  OF  THE  FACTORIZED  H»S>*> 


-<w 


— 


1 


*  ; 

■M 


CALL  PRCVEC  (C  A,  Nl 
C 

C  CONVERT  THE  FIRST  ORDER  PARTIAL  F  f ACT  IONS  TO  7  CCMATn 
C 

NR  T  HO  =  IHT  HD  +1 

GO  TO  C2C30.300C.400CI.  NHTHO 
C 

C  IMPULSE  INVARIANT 

c 

2000  002  3091  =  1. N 

CA  (  II  =  C4  ( II  »OEL  TA 
20 C  9  CR  <  II  =  CEXC  (Ck  (I  »»0FLTAJ 
GO  TO  45  0  3 
C 

C  PULSE  INVARIANT 

C 

3000  003309  1=1, N 

C0N1=C£XP <CR(  I*  ’DELTA! 

CACI!=CACI!*(CON1-1.0I/CR<I) 

30  C  9  CR  ( I ) =C0N1 
GO  TO  4500 
C 

C  TRAPE2CI0AL  INVARIANT 

C 

40  00  I  CM  ECK  =  2 

004009  1=1, N 
C0N1=CEXP(CR(II*0£LTA> 

C0N2=CA(I>/(CP(I>’CMII*CELTA*C0N1I 
CONT=CONT*CON  2*  C(1.3-CR< I>*OFLTAI  *C0N1-1.  C> 

CACI»  =  CONO*<1 .0  -CONI  >*<1.0 -COM) 

4009  C R  <  1 1  =  CON  1 
GO  TO  4500 

c 

c  CONSTRUCT  THE  Z  OOMAIN  DENOMINATOR 

C  ANO  NUMERATOR  POLYNOMIALS 

C 

4500  CONTINUE 

CALI  PCST7<CP,C8,0,NI 
00  4509  1=1, N 
45  0  9  CA  A  (I)  =3.  COO? 

00  4519  K=1,N 

CALL  PCST7(CR»CA1,K,NI 

00  4519  J=1,N 

4519  CAA «J) =CAA ( J) *CA1 (J)*CACK) 

CAA  ( NP 11=3,2 

C  ****•• 

IF (IMTHQ.NE.l >  GO  TO  4521 

DC  4523  1=1, N 

II=NP1-I 

45  20  CAA(II»1I=CAA (III 
CAA  (H -0. 0 
4521  CONTINUE 
C 

C  AOJUST  FOR  DIRECT  TRANSMISSION 

C 

C  DTXC=(fl. 3,0,2 I 
C  IF(NN.N£.n|  CONT=  DTXC 

C  CAA(NP1)=C0NT»CB«NP1I 

00  4529  1=1, N 

4529  CAA  (II =CAA(I) ♦C0NT»C8(I» 

C 

C  SHIFT  NUMERATOR  TO  COMPLETE  PULSE  INVARIANT  TRANSFORM 

C  WHEN  NUMERATOR  HAS  LOHEK  orcer  than  denominator 

c 
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CONTINUE 


C 

c  - 

c 

c 

C  PRINT  THE  TRANSFORMED  CCEFFICIENTS 
C 

50  CC  CONTINUE 

c 

WRITE (6,55101 

95 1C  FORMAT (*  POLES  IN  T«E  Z  DOMAIN*) 

CALL  PRCVECCCR.N! 

HR  I  TE  ( 6 , 55  20  FAC 

9520  F CRMAT  (*  GAIN  FACTOR  USEC  =*,E14.7,//) 

WRITE(6.953C) 

9530  F ORMAT  <*  ZEROS  IN  THE  Z  DOMAIN*! 

CALL  PRC  Vr C  (C  A,  NN ) 

WRITE(6,9*4C) 

954  0  F CRMAT ( *  Z-OOMAIN  DENOMINATOR*! 

CALL  PRCVEC(C8,NP1I 
WRITEI6,955fl> 

9550  FORMAT  < *  *-OONAIN  NUMERATOR*! 

CALL  PfiCVEC<CAA,NPll 
00  5015  I=1,NP1 
8(11*00(1) 

5515  A ( 1 1  =  C  A  A ( T ) 

GO  TO  113 
C 

5029  WBITEI6. 95631 

956  C  FCRMAT  (/, IX,* NUMERATOR  ORCEF  LESS  THAN  ZERO*,//! 

5999  STOP 
EN,1 

SUBROUTINE  PCSTZ(C»R2»K, M 
C 

C  FCSTZ  CONSTRUCTS  A  Z-COMAIN  POLYNOMIAL  COEFFICIENT  ARRAY 

C  from  AN  ARRAY  OF  ITS  ROOTS. 

C 

OIMENSION  CU1.R2C1! 

COMPLEX  C,R2 
NP 1=N* 1 
001I=Z,NP1 
1  R2(I1=0.0 

P2(1I=1.C 
0031=1, N 
IF  (  I-KI  6,  5,6 
6  004  JJ= 1,1 

J=I-JJ*1 

4  R2(J*1I*P’(J*11-C(I)*PZ(J) 

3  CONTINUE 

RETURN 
ENO 

SUBROUTINE  PRCVEC  ( A,  Nl 
C 

C  PRCVEC  PRINTS  A  COMPLEX  VECTOR 

C  C  A  COMPLEX  NUMBER  OF  THE  FORM  A*  J8  IS  FRINTED  (  A,  8  J) 

C 

DIMENSION  All! 

COMPLEX  A 

IF(N.EQ.O)  GO  TO  100 
WRITE (6,920! 

WRI TE 16,910  I (A( II, 1=1, Nl 
91 C  F0RMAT(1X,IH(,F22.15,1H,,F22.15,3h  Jl! 

100  WRITE (6, 920) 

92 C  FORMAT (21/)) 
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c 

c 

c 

1 

31 


c 

c 

c 


1C 


2 


4 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


return 

END 

SUBROUTINE  PRVEC<A,N> 

THIS  SUBROUTINE  OUTPUTS  A  SINGLE  DIMENSIONED  ARRAY 
DIMENSION  A(1 » 

WRITE<6,1>  <A< I) ,1=1. N) 

FCRMAT<lX,10Fl3.5t 
FORMAT (/> 

RE  TURN 
END 

SUBROUTINE  FOLCON  CC.  R2  , K  , N  > 

A  POLVINO'dAl  CONSTRUCTION  PROGRAM  NEFQEO  F CR  Z 10  S 

DIMENSION  C  U  >  ,R2  <11 
Complex  c,r2,comp 
REAL  OC<?l 

EQUIVALENT  CCOMP.OCl 

NP1=NM 

00101-2, N°1 

R2<II=3.C 

R2  <lt  =  l.G 

00**  I  =  1  ,  N 

COMP  =  C  C I > 

IF < T.EQ.K.OP.  <OC<l).EQ.O,C.ANn.Or<2l .EO.C .0  l»GO  TO  „ 

302 JJ=1,I 
J=I-JJ»1 

<?2  <  JM)=R?  <  J*  11  »C  <I>*R2I  J) 

R2<l*  =  P2<l»*C<n 

continue 

return 

ENO 

SUBROUTINE  POLRTt  XCOF.COF, M,R OOTR , ROOTI . I ER » 

COM°U  TE  S  THE  REAL  AND  COMFLEX  ROOTS  OF  A  REAL  POLYnO‘Hl 
DESCRIPTION  OF  PARAMETERS 

XCOF  -VECTOR  OF  M*1  COEFFICIENTS  OF  THE  pCLYn'-mTAL 
OROEREO  FROM  SMALLEST  TO  LARGEST  PO HE R 
COF  -WORKING  VECTOR  CF  LENGTH  Mfl 
M  -ORDER  OF  POLYNOMIAL 

ROOTR-RESULTANT  VECTOR  OF  LENGTH  K  CONTAINING  RCA t.  FOOTS 
OF  THE  POLYNCMIAL 

ROOT I-RESULTANT  VECTOR  CE  LENGTH  M  CONTAINING  THE 

CORRE SPONCINC  IMAGINARY  ROOTS  OF  THE  POLYNOMIAL 
IER  -ERROR  CODE  WFEPE 
I SP=C  NO  ERR  CR 
IER=1  M  LESS  THAN  ONE 
I ER=?  M  GREATER  THAN  36 

IE®  =  3  UNABLE  TO  CETERKINE  ROOT  WITF  5K  INTE RATIONS 
ON  5  STARTING  VALUES 
IER=<*  HIGH  ORDER  COEFFICIENT  IS  ZERO 

DIMENSION  X  CO  E  < 1 1 , COF <11 ,R00TR<1I  ,ROOTI Cl > 


LIMITED  TO  36TH  ORCEP  POLYNCMIAL  OR  LESS. 

FLOATING  POINT  OVERFLOW  MAY  OCCUR  F CR  HIGH  CRDER 
POLYNOMIALS  BUT  WILL  NOT  AFFECT  THE  ACCURACY  OF  THE  PfSULT 

METHOO 
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NEWTON- RAPHSOP  ITERATIVE'  YE CMNI QUE.  THE  FINAL 
ON  EACH  ROOT  ARE  PERFORMED  USING  THE  ORIGINAL  POLVHP-IAL 
RATHER  THAN  THE  REDUCED  POLYNOMIAL  TO  AVCIC  ACCJMULtT-0 
ERPORS  IN  THE  REDUCED  POLYNOMIAL. 


IFIT*C 
N=M 
IER  =  0 

IF«XC0F«N+1»I IP. 25, IS 
10  IF(M  15,15,32 

SET  ERROR  CODE  TO  1 


15  IEP=1 

20  IF  (  IE  R  >  2  3 1  .  2C  1 .  2  0  C 
} C  WRI TE  <  6, 2S3 ) I E® 

33  FORMAT  (IX, ’ERROR  CALLED  FROM  POL®  T ,  IF.o  =  •  .  1 31 
31  RETURN 

SET  ERROR  CODE  TO  4 


25  IER*** 

GO  TO  20 

SET  ERROR  CODE  TO  2 


30  IE R  =  2 
GO  TO  20 

32  IF ( n-361  35.35,3- 
35  NX  =  N 
N  X  X  *N  + 1 
N2=  1 

KJ1  =  N’l 
OO  40  L=1,KJ1 
MTsKJl-L’l 
40  COF  (MT l=XTOF(L> 

SET  INITIAL  VALUES 

45  XO= , 0C  5  K 10  1 
YO*0.310CC101 

ZERO  INITIAL  VALUE  COUNTER 


IN*fi 
50  X*XC 

INCREMENT  INITIAL  VALLES  ANC  COUNTER 

XO*-1J.O*YO 
YO*-10 ,0»X 

SET  X  ANC  Y  TO  CURRENT  VALUE 

x*xc 

Y=  Y0 
I N* IN’  1 
GO  TO  59 
55  IF  I T  =  1 
XPR  =  X 
YPR  =  Y 


EVALUATE  POLYNOMIAL  ANC  DERIVATIVES 


64 


o  o  o  o  o  o  ^ 


59  icr=o 

69  U  X  =  3  •  i> 

UY=|].J 

V  =3.0 
vT=o.a 
X  T  =  1. 0 

U  =  CCF  <N,1< 

I ^  <  U *  65,130,  65 
65  00  70  1  =  1. N 
L  =N-I*1 
TEMP=C0r<i.» 

XT2=X*  XT- Y* YT 

Y  T 2  =  X*  »T*v»XT 
U=UyTE»*o»TT2 
V=VyT£mP*yT2 
FI  =  I 

l)X  =  UX*PI»XT*T  EH  P 
UY  =  6Y-FI»-'T*T£HP 
XT  =  XT2 

70  YT  =  YT2 
SUNSQ=UX*UX+UY*UY 
XF(SUHSO)  75,110.75 

75  0X= <V*UV-U*UK 1/SUHSO 

x=  x*ox 

OY=-<U*UY*V*UXI /SU*SQ 
Y=  Y ♦ 0 x 

9  IF  <  ASS  <1Y»  ♦  A«  S<OX  l-l.CE-ltl  1C  C.  9G  ,80 

STEP  ITERATION  COUNTER 

93  ICTsICT*l 

IFJICT-50:i  60.95.85 
85  IF«IFIT*li?.90.13C 
90  IF<rN-5»  JO. 95. 95 

SET  ERROR  COCE  TO  3 


95  I E  R  =  3 
C.0  TG  2J 

103  DO  105  L=1.NXX 
HT=KJl-l+l 
T E H P=  XCQF  (HTI 
XCOP(HT>  =  rOFH> 

1C5  COF (L> =TE yP 
I T  E  HP=  N 
N=NX 
NX=IT£MP 

I F  ( IFI T I  1  20, 55  ,12C 
110  IF ( IFIT1  115,50,115 
115  X  =  X  PR 
Y=YPP 
120  I F I T=0 

122  IF  ( ABS  « Y »  -1.0  E- 8*  ABS  (XI)  135,1  25,125 
125  ALPPA=X*X 

SUHSO=Y*X»Y*Y 
N  =  N-2 
GO  TO  14C 
130  X  =  3.0 
N  X  =  NX- 1 
N  X  X  =  NX  X- 1 
135  Y  =  0  .0 

SUHSO=  0.0 
ALPHA=X 
N  =  N-1 


_ _ _ ,■  y 
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143  COF(2>«COF|2>*ALPPA*COFU> 

145  00  153  L=2,N 

15C  COFa*l*  =  roFJL»l>4ALPHA»CCFai  -SUHSQ»C0F(l-1> 
155  »00TI<N2»=Y 
S00TRIN2)=X 
N?*  N2* 1 

IF<SU*SQ>  163. 165,160 
160  Y  =  -Y 

SU43Q= 0 . C 
GO  TC  155 

165  IF(M  23,20,45 
ENQ 
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• •••* »»»«»»»»•»»•*,***«»•♦»»»» 

PROGRAM  "GOUANT“ 

INPULSE-OE'PCNSE  MODELING 
BY  P£  NC IL-CF -FUNCTIONS  METHOD 
DEC  1979  «FGR  RACGI 

PROGPAM  "GQUANT”  USES  CHARACTERISTICS  OF  OUANTI7ATION  £R=C~ 
IN  PENCIL-CF -FUNCTIONS  hethOO  TO  PROCURE 
IMPROtfEC  TRANSFER  FUNCTION. 

MOOELS  IMPULSE-RESPONSE  OF  CHA KN EL /N -I NO R< 
CAN  BE  USED  IN  SINULAT ION  MCDE  <eiN.  0°  Cr 
OR  ON  EXPERIMENTALLY  RECORDED  CESPONSES. 


OIMENSION  F  <E  00  l.U  <50:  1 ,  LU  <5CC  I  ,X  <51:  ,81  ,  G<  P.  «l  ,AM<8,8| 
01 MENSI ON  GN<4,4I,G-ST<9,8I ,GCUM<8 ,81 ,E  <8 , 6  I  »  EN  <  H  ,  *  | 
DIMENSION  V<l6l ,VV<16I ,AMF  <»l ,SR<PI ,SI(8I ,SPH<8 I 
OIMENSION  TITLE<70),I8UF<512I 
DOUBLE  PRrC IS  ION  OT, AC, 30, ERROR 

COMMON  /04j/ISPN,XMS9,0ELTA,SIG2,CT,0I,BIAS .IBI AS 
COMMON  /0A1/FB4R, E8AP, F£ SUM, FESUM 
COMMON  /IO/IR.ILT, IPR, IRCUNO, IPLT 
RERIN05 
MAXPL=50C 
M4X  =  8 
MAX2=2*“A* 

IR  =  5 
ILT =6 
ISKIP=0 
NSTST=2 

CALL  VEQUAT (MftXPL ,U,F.O. 1C» 

CALL  YEOUAT  IM  AX  PL  ,  UU ,F ,  0 , 1 1» 

CALL  WE0U1T<mAX2,V,VV.:.C» 

MRITE<ILT,21 

REA0<IR,8I  <TITLE<II,  1  =  1.70 
MRITE< ILT, 181  ITITLE<II ,I  =  1,7CI 
REAO  (I  R,81  <  TITLEdl,  1  =  1, 70 
REAC<IR.8i <TITL£<ll«I=l,7CI 

REAOeiR.LtNPT . IRAD.NO IG.N, ISIM,NCCMP, IPLT  , NNPT, 
♦XNSB.OT.BIAS 
NP1=N»1 
NP2*NPl*l 
NP3=N*3 
NPNP2=N+N<2 
NPNP1=N*N*1 
IFINNPT.EO. 0) NNPT  =  NPT 
IF  <CT. £0. j. 01  0T=1 . 0 

IF  <  ISIM.EQ.  Cl  REAO  (TR,  18C  )  <F  <K I  , K=  1 ,  NNPTI 
IF |  ISIM.EQ. C»GO  TO  61 

IF  (ISIM.EQ.  It  RFAO  <IR,19C1  f  V  <1 »  ,  1=  1 ,  NPNP2I 
TF(  ISIH.ET.il  CALL  RESPON  (F,U,N,V,VV » NNPT I 
IP ( ISIM.EQ. II  GO  TO  61 
00  6C  I=1.NC0MP 

PE  AD < I R, 5) AMP  <11,  SR  (II .  El  III.SPMII 

60  MR  I  TEC  ILT,  111  I,  AMP  <11  ,  SR  <  1 1  .  S  I  <  1 1 ,  SPH  <  1 1 
CALL  SIGNAL<F,NNPT,AMP,SK,SI, SPM.CT.NCOMPI 

61  CONTINUE 

IF  UPLT.GE. 21  CALL  PL OTS < IBUF, 5 1 2. °> 

1111  REACCIR.8I <TI TLE(II , 1=1, 7£ I 
IF<EOF<IRI.NF.OIGO  TO  99 € 
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WRITECILT.31 

WRITE  <  ILT. 181  (TITLECII ,1  =  1,70  I 

RE*G<IR»4)IPR»IRE7,ISPK,IFIX,KFIX,I8IAS,IYV,IZ7,CI 

IROUNO=d 

C  ROUNDOFF  3PTI ON 

C 

41 0  CONTINUE 

IFCIROUNO.NE. 0) CALL  Ou  AN  T  Z  CF ,  X,  NO  I G ,  IF  AC. NPT,»A  *PL I 
88  CONTINUE 

IF(IPOUnO.NE. 01  GO  TO  99 
00  30  K=1,NPT 
31  X OC.ll =F<«J nBIAS 

99  CONTINUE 

IF (NP1  .GT .11CALL  INGPATCX.NPT  ,NP1 , m A XFL  , - 1) 

C 

C  CONPUTE  G®  AM  MATRIX 

C 

NPP=NPt 

IF ( I0IA3.NE.C )NPP=NP2 
00  44  1=1,  NPP 
00  44  J=1,NPP 
AO=C.O 

IFdSPN.EO.O.  AND.  jrcund.fc.  :)  GO  TO  -3 
DO  42  *=NSTRT,NFT 

42  AD  =  A0+XCK,I1*  XC  X,  Jl 
GN«I,JI=A  l»OT 
GOUMII, J»=GN<I, J> 

43  CONTINUE 

IF< IROUNO.EQ. Cl GCI,J)  =  GN  CI.JI 

44  CONTINUE 

rFdSPN.NE.C.OR.IRCLNC.NE.?! 

1CALL  GKRO'.T  <gn,  E.CET,  V,NFF,NPP,IA»,1) 

IF ( IPO  UNO. EO. Cl  WRITE  CILT, 171) DET 
IF  (IRCUNO.EO.  II  W®  TTE  CILT  ,1721  OET 
IF( IPR.GE.1ICALL  P»TMAT<GN,NPF,NPF,MAX,-1  I 
C  wRITEdLT.il 

IPO  =  IROUN  i 

IF (IRCUND.EQ. Cl I®CUN0=I®CLN0»1 
IF C I®0.E3. 9 .ANO.TSFN.NE. -1) GO  TO  410 
IF ( IFIX.EO.-l IGC  TO  203 
C 
C 

C  ESTIMATE  OF  ••  G 

C 

156  CALL  8UIL0Z (AM,V.NPl,NPT,MAX,NFIX) 

C  - N P 1  REPLACED  BY  HP ®  NEXT  3  CARDS-” 

CALL  FIXIGOUM ,AX, GEST.E, V, NFP,NPP,SIG2,MAX, IF  IX 1 
IF  <  IFIX.EO.il  WRITE  (ILT, 4821  SIGZ 
CAUL  GXROCT CGEST, E, OET, W, NPP, NPP, MAX , II 
WRITE! ILT, 162I0ET 

IF C IPP.GE.1IC  ALL  PPTM AT CGEST, NP 1, NP 1,M AX , 01 
00  154  I=’,NP1 
00  154  J= 1 , NP 1 
154  GOUM(I, JI=GEST< I, J1 

NF I X=NF I X- 1 

IF  CNFIX.GE.il  GO  TO  156 

ls«dP=i 

IROUNO=0 

c 

C  CALCULATE  ERPCR  MATRIX 
C 

150  IFCIPR.LE.2.0R. ISKIF.E3.0  ICO  TO  151 
00  32  1=1, NP1 
00  32  J=  1 , NP1 
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E(I,J)=G(I, J> -G£ST(I, J> 

32  EN(I,J)=G(I,J)-GN(I,JI 

MR  I TE  ( ILT  ,  1 61  > 

CALL  PRTNAT  (E  , NF1 . NP1 ,NA X  ,-l> 

MR  I  T£  (  ILT ,  1  63  > 

CALL  PRTMA  T  (EN«NP1*NP1,MAX,-1) 

151  CONTINUE 
C 

C  DETERMINE  NUMERATOR 

C 

203  CONTINUE 

IF(ISPN.EC.C>  GO  TO  998 

CALL  VEOUAT (NPl.V (NP2 ) ,V V, 0,10) 

CALL  RESP1N(X (1.1 1 ,U,  ,V,VV,NPT) 

CALL  ING3AT (X , NPT , N°1 -IR i  * , PA XPL.2I 
C  CHANGES  PAOE  HEREON  FOR  6  €  2 » - C 

L=N-IR£N 

IF(  IBIAS.NE  .0  1L=N-IRFP«-1 
L  PI  =L  ♦  1 
L  P  2  =  L*  2 

IF (IBIAS.NE.c (CALL  VEOUAT  ( NPT , X (1 , LP1 » , U, f . 11 > 

CALL  VEOUd  (NPT,X  (1.LP2I  ,F,0,  1> 

CALL  VEOUAT  (NPT,X (1,LP2) , PI  AS,  0,3 > 

DO  21b  1=1, L 
DO  216  J=  1 ,  LP 1 
G(I,J»=C. . 

DO  215  <=l,Nt>T 

215  G  (I,J>=G(I,J>  4X(K,I*ll*XO<,J*ll 

216  G(I,J»=G(I,J»»OT 

C20 5  CALL  PRTMATiG,L.LPl,MAX,2C5> 

CALL  G<ROCT(G,E,OET,VV,L,L,NAX,C> 

C207  CALL  RRT-AT<  E,  L,  L ,M A  X,  2  0  7> 

CALL  VEOUAT <NP1, VV.AKP.C  ,C> 

DO  219  1=1, L 
DO  219  J=i  ,  L  •• 

219  VV (I»=VV(I»+£ (I,J»*GCJ,LP1»/D£T 
fmean=:.o 

IF ( I  PI AS.NE.C  »FMEAN=VV (L  > 

V(NF2 1=0,0 

CALL  VEOUAT <N,V(NP3» ,VV. C.l) 

211  MR  I TE ( ILT  * 3C?  > 

MRIT£(ILT.21C I (V(I»,I=1,NP11 
M  R I T£ ( ILT , 2 10 » (V(I»,I  =  NR2,NFNP2» 

IF ( I9IAS.NE.0  >MRITE(ILT, 3C5»f»EAN 
CALL  RESPON  (X  (1 , 2 1  ,U , N ,  V ,  V V,NNPT> 

ERRGR=a.C 
FFSUN=0,0 
DO  213  K=  1 ,  NN FT 
FFSUH=FFSUN*F  «M)*FfKJ 
X (K,3I =F (K» -X (K,2  » 

X  (X.l)  =F(M 

213  ERROR=ERROR*X (K,3)*X(K,3> 

FFSUM=FFSUM*DT 

ERKCP=FRROR*OT 

RATIO=£RROR/FFSUM 

MR  I  TE(  ILT,  304  »E  RR  CR,  FFSUP.  RAT  IO 

IF (IPR.G£.2»MPITE<ILT,11C» (FCIO  ,K*1,NNPT» 

MRITE(ILT,1> 

IF(IPR.GE.2IMRITE (ILT, 1181 I XI X,  21  ,  K=i  ,  NNPTI 
WRITE (ILT, II 

IF( IPP,GE.2»HRITE(ILT,110I (X(X,3» ,K=l,NNPT» 

C 

TO  =  C  .  0 

IF(  IPLT.GE.  2»  CALL  PLOP  (NNPT  ,2  ,  X,M  A  XPL  ,T0  ,  CT  ,  1HV  ,  1HT  ,  I  BUF> 


format  statements 


c 
c 

5  FCRMAT(5F1C .3) 

4  FOR  MAT  (I*«,6I2,l4,4FlO.Ci! 

6  FORMAT (70  A1 1 

IS  FORMAT (2X.70A11 

903  FORMAT  (10  (5*.  F5.C  )  > 

93  4  FORMAT  (1C  (5*.  15!  1 

11  FORMAT  ( 2X , 1 2,  ’  Am  P  =  »  ,  F  8  «  2  ,  *  S=*,F10.4.*  ♦  J*,F1C.*-, 

1*  PWASr  =  ».FlC.<,| 

14C  FORMAT (1CX, 8HG  MATRIX! 

16 0  FORMAT (10 X. SHM  MATRIX* 

161  FORMAT  (1C  ,  12  HO.  TRUE  -  G'ST> 

163  FCR“AT(10x,1hHGTdUE  -  GNC1SY) 

162  FORMAT  (10X, 11 HGEST  “ATR1X,*  (DE I  =  * ,G1 3 . b  , * ) • I 

210  FORMAT  (7X, 5 (2X.G1 3,61 ) 

110  FORMAT  (2G  <1X.  F5.2  »  I 

C11C  FORMAT  (2X.10(  2», FIS. 5» » 

178  FCRMAT(10x,14HNOISm  X  rATPIXI 

179  FORMAT (li»,8HX  MATRIX! 

171  FORMAT  (1CX,  16  l-TRUE  GRAM  MATRIX,*  (3£T  =  *  ,C- 13  .6  .*  I*  « 

172  FORMAT (1QX,17HNCI  1 Y  GRAM  MATRIX,*  t  OET*  *,  G1  3 . 6,  *  I  *  > 

18  0  FORMAT  (5F10.C 1 

30  3  FORMAT  (2X,  *EST  NL'M/O-MOM  VECTOR*  1 

33  4  FORMAT  (/,  X,*  fcRPOR=*,Gl3,6,*cFSUM  =  *,G13,6**AATIC-*,Gl  i.r>,  /! 

335  FORMAT (2X. 'ESTIMATED  MEAN* * . G 1 3 . 5 > 

432  FORMAT  (*  ESTIMATE  NOISE  VAR**,Gi?.5l 

1  FORMAT (/! 

2  FORMAT  (lm  I 

3  FORMAT (////) 

C 

C  1 5  I M=  0  FO1-'  MODELING  ACTUAL  »riR OUSE  DATA 

C  1  OR  2  FOR  SIMULATION  (ll“(7!l  (2tSUMS  CF  EXP  AnC  OSCl 

C  NO  I G*0 IN  PITS ( I NC  SN-8IT1,  OR  DECIMAL  MANTISSA.  ROUNDOFF  In  a 

C  IRAC  =  2  Fqy  dlNARx,  1C  FOR  DEDI*.-L 

C  NPl=INTEGrATED  FUNCTIONS,  THE  FI^ST  IS  DATA 

C  I  PR  =  C  FOR  MINIMAL  PRINTING,  OTHERWISE  1  OF  2 

C  ISPN  =  C  TF  ANALYSIS  OF  R.CFF  FRFOP  SIGNAL  CNLY.l  FOR  Pf.JN.T"'  S 

C  -1  IF  ANALYSIS  OF  TRUE  (LNPOUNDED*  SIGNAL  ONLY 

C  NCOHP=  COMPONENTS  (A*f  XP(JR  T!  *  SIN(SI  T !.!  TYRE 

C  INT=1  (OR  31  FOR  FORWARD  INTEGRATION,  -1  FOR  REVISE 

C  OT=SAMPLlNG  INT  XMse=W£IGNT  OF  MS*  <*(N-11! 

C  NPT  =  OATA  PTS. ,  NNRT  =  POINTS  ON  PLOT,  N=0  EG  RE  t  Cc  HCO^L 

C  IF  I X  =  * 1  FOR  NO  CORRECTION,  1 

C  NF I X=1  FOR  07  CORR CC • ,  2  FOR  OZ  Q7  CCRREO 

C  3  FOR  81  AS  ANC  Q 7  07  “OPRFCTION 

C  I F I  x=  0  IF  GEST=GN-Am,  i  IF  NCISF  VAR  TO  -PE  ESTIMATED 

C  -1  IF  NC  CORRECTION  IS  TO  8E  APPLIED  (ISFN  MLST  eE  11 

C 

GO  TO  1111 
9*8  CONTINUE 

CALL  PL0T(C..C.,999! 

STOP 

ENO 

SU9RQUTINR  SI  GNAUF.NPT,  AMP.SR.SI  ,SPM,  DT,  NC  CM  R1 

C  - - - 

DIMENSION  F(1 I , AMP(i l,SP(ll ,SI(1I ,SPH(1I 
COMMON  /IO/IR,ILT,IPR,IRCUNO 
OOU8LE  PRECISION  A,8,C,X 
00  12  K=1,NPT 

12  F (KI  =  J  ,0 

00  20  1*1 , NCO  MR 

A=SR(I!*OT 

9=31(11*07 


70 


C=SPH( I) 

00  15  KK=1,NPT 
K=KK-1 
X=AMPI I) 

IF (A.NE.S .0)X=X*OEXP  <A*K> 

IF  <  6.NE.0. B)X=X»OSIK(€»K+C> 

15  F(«K)=x+F<KK) 

20  CONTINUE 

c 

IF<IPR.LT.2)G0  TO  3C 
WRITEULT.9) 

WRIT£<ILT,6  > <F<  K) ,K=1,N°T) 

WR  I  TE  ( ILT  ,  1 1 
30  CONTINUE 

1  F  CR  PAT  (/ ) 

6  FCRPAT  (2!M1X,  F5.2  H 

9  FCRMATUOX,*  F  SIGNAL  * ) 

RETURN 

ENO 

SUBROUTINE  CUAN  TZ  <F,  X,  NOIC-.  IP  AD.NFT  ,NDIN> 

C  - - - - - . . 

C  PERFORMS  BINARY  OR  DECIMAL  QUANTIZATION 

DIMENSION  F  II  l,X  (NOIM.  1) 

DOUBLE  PRECISION  ct.ac.bc 

common  /OAG/ISPN, XMSB, DELTA, SIG2.CT .01 , 31  AS , 1  El  AS 
COMMON  /DAI  ABAR.EBAP.FESUM.EESUM 
COMMON  /IO/IR.ILT.IPR.IROUNr 

c 

C9 

FflAR-g  * 

EBAS=0 . 

FESUN=C. 

EE  SUM=  0 • 

C 

C  BINARY  QUANTIZATION 

C  WO»0=SN  BIT.MSB,.  ...LSB  MAX  NEG=-2*XFSB 

IF (IRmD.NE.2) GO  TO  551 
IF ( XMSB.EQ.G. C>  XMSB=5.0 
N0IG1=NQI  -  1 

OEL  TA= <2.G*XM  SB)/ C2.C**NCIGl> 

0EL=0ELTA/2.C 
SIG2=OE».Ti*OELTA/12.  C 
W°ITE<ILT.**B9)0ELTA.5IG2 
DO  81  K=1 .  NPT 
XLEV=i  .0 
SN=-1.0 
XX=F(KI*8TAS 
IFIXX.GE.j.C) SN= 1 . 0 
XX=SN*XX 
00=2.0*  XMSB 
X0=2.C*XHSB-D ELTA 
00  82  1=  1.  N  01  G 
00=00/2.0 
OIF=ABS(XX-XlEV) 

IFICIF.LE.Oa  )XQ=XLEV 
X  TEN=XLE  V-00 

IF  I XX.GE.XLEV)XTEM=XLtV*CO 
XLE  Y=XTFM 

C  WRITEIILT .2  1C )XX. DC.OIF.XLEV.XO 

82  CONTINUE 

IFIXX.LT. DEL)  XQ  =  0  •  G 
XfK,l)=SN*XQ 
81  CONTINUE 

GO  TO  711 
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c 

Cm2 

CAfc  OECIMAL  O'lANTIZATICN 
551  CONTINUE 

JF<IRAO.N£.10 1G0  TC  711 
AAA=lu .G**NCZG 
0£LTA=1. 0/AAA 
SIG2=0ELTA»0£LTA/12.  f 
WRITE (ILT  .  A89 I0ElTA,S1G2 
00  91  K=1,NPT 
X(K,1»=G. ' 

FB*  F  CKI *61  AS 
XTEM=AES (cBI 
SN=1.G 

IF (FB.LT.:. Cl  SN=-1.C 
XT£M=XTSM*  AAA 
XT£M*XT£M+G.5 
I X=  XTEM 
XTEM=IX 
XTEM=XT£N/AAA 
X (K,1)=3N»XT£N 
91  CONTINUE 

C 

C  SSO  VALUE; 

C 

711  CONTINUE 

00  211  X=1.NPT 
F0=F(X) .BIAS 
X  CK.2I=X  CK. II -Fp 
FBAR  =  FBA<J*FB 
EBAR=EEAR*X (K,2> 

E£  SUM=  EE SUM  .X CK.2I *X  <x, 21 
FESUM=F£$UM.F  B»X(K,21 
IFCISPN.EO.CI X(K. 1»=X(K,2» 

211  CONTINUE 

eesum=eesum»ot 

FESUN=2.C*FESUM*0T 
FBAP=P0AP/NPT 
E  P  A  8  =  E  8AR/  NPT 

WP ITE  < lLT,Na2  >FBAP,E8AR, FESUM .EESUK 
IFCIPR.LE.2I  GO  TO  *.11 
WRITEdLT.  8) 

MPITECtLT.llO  MX(K,H,k=1,NOT> 

IF ( ISPN.EQ. U GO  TO  All 
WRITECILT.  1*1 

MPITEC ILT.115 I CXCK.2I ,K=1,NFT> 

WPITEC  IlT.1I 
All  CONTINUE 
999  CONTINUE 
C  FOR  PAT  STATEMENTS 

C 

8  FORMAT <1G*. •ROLNOEC  F.glAS  SIGNiL»l 

18  FORMAT  C 1 0 X . OUNOOFF  EPROP  E*  ) 

21  C  FORMAT  <2X  ,5  (2  x,  Gil  .A)  ) 

11C  FORMAT (20(1X, F5.2  II 
115  FORMAT (IX, 20( 1X.F5.3I I 

178  FORMAT  C10X,  14  l-NCISV  X  MATEIXI 

179  FORMAT C1CX.8MX  MATRIX  I 

A82  FCRMAT<2X,5MFBAP=,F11.A.EM  EB AR*. El 1 . A, 5P  FE2  =  , El  1 . *. LH  tr*,tll.-l 
A9 9  FORMAT (2X,6«3ELTA=,Fi;. J, 5m  SIG=.E12.AI 
1  F  CRMAT (/I 

RETURN 
ENO 

SUBROUTINE  COPUPT (F, X , NO IG, IP AO.NPT , NOI “I 
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c  - 

c  aoos  Ncrs- 

0IMEN3I0N  F (1 ) ,X(NCIM,1) 

OOU6LE  PRECISION  OT.AO.PC 

COMMON  /OAO/I  S°N,  *MSB,DELTA,SIG2, OT  ,QI,8I  AS  ,IBI  AS 
COMMON  /0« 1 /F BAR,EeAP,FESUM»t ESU« 

COMMON  /IO/IP .ILT.IPR.IRCtNn 
C 
C9 

FBAP=3 . 

FBAR--0. 

FESLM=a. 

EESUM=C. 

c 

NOIGl=NOIG-l 

0ELTA=(2.;*)(VS9>/<2.C**NCIG1) 

SIG2=OELTA*OELTA/12.C 

WRITE(IlT,489>0ELTA,SIG2 

c 

c 

15=295816= 

IS2=397665 
SIGMA=SQR* (SIG?) 

CALL  NRML (NP7, 1,1,0, , SIGMA, IS  ,IS2,X (1,21, C» 

DO  26  K=  1  ,  N  FT 

26  X(K.l>sF(K)«SIAS*X<K.2> 

C 

DO  211  K= 1 , NP7 
FB=FCKM8IAS 
F  B A  R=FB AR  *F  B 
EBAE=£9AP*X  (K,2) 

EESU«=ESSUM*X  (K,?)*XCK,2> 

FESUM=FESU,1*F  P»X(K,2> 

IF  (  ISPN.E9.  lit  X(K,  1I  =  XIK,2» 

211  CONTINUE 

EESUM=EESUM*OT 

FFSUM=2.G*FESUM»OT 

FBAR=FBAR/NPT 

EBARsEeAR/NPT 

WR ITE ( ILT ,4  82  >FBi S  ,E e A R, FE SUM . EES UM 
IF(IPR.lE.2»  GO  TO  911 
NPITEIILT.8) 

NPITEI ILT,  110  KXIK.1 >  ,  K=  l.NPT  I 
IF ( ISPN.EO. C) GO  TO  911 
WRITE  ( ILT.  181 

WRITE (ILT, 115  I (X«X,2) ,K=1,NP*» 

WRITEdLT.il 
911  CONTINUE 

999  CONTINUE 

C  FORMAT  STATEMENTS 

C 

A  FORMAT <lOX,16MROUNCEO  F  SIGNAL) 

18  FORMAT IlOX, 16WROUNCOFF  EFROR  E) 

210  FORMAT (2X.5I2X, Gil. 91) 

11C  FORMAT  <23UX,F5. 2  I) 

115  FORMAT (1X,2C( 1X.FE.3)) 

982  FORMAT  (2X,5)FBAR=,E11,9,6H  EBAR  =  ,E  1 1 .9, 5H  F  E2=,  El  l .  i,,  9«  EE=,E11.9) 

989  FORMAT  <2X,6H0ELTA=,FU  .3, 5M  SIG=,E12.4> 

1  FORMATI/I 

RETURN 
ENO 

SUB  ROUT  IN  -r  ingrat  ix.npt,  NP1.no  I  M.INT) 

C  - - - - 

DIMENSION  XINGIM.l) 
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OOUBLE  PRECISION  tT.SC.9C 

COMMON  /D-  O/ISRM.  XMSR.OELTt.SIC2.CT.OI.SI  AS  .161  AS 
COMMON  /IO/IP  *ILT . IPP, I R  OUN  C 
C  GENERATE  INTEGRATEC  SIGNALS  FROM  CAT6  IN  XtK.ll 

C  INT  =  1  OR  2  FOR  FORWARD  INT.,  -1  FCR  REVERSE 

C  I  NT  =2  FOR  UNIT  DELAVS  <X  <K,  !♦  1»  =X  IK- 1 , 1 )) 

N=NF1-1 

NP2=NP1*1 

I0PT=INT»2 

GO  T0<51,ll,ll,9i»,TGPT 
C  FORWARD  INTEGRATION' 

11  CONTINUE 

DO  AO  J=1.N 
JJ= J+l 

*(1,JJ>=X<1,1> 

00  AO  K=2.N°T 
Kl=<-1 

XIK, JJ»=X  <K1. JJ>  +X  CK,  J» 

AO  CONTINUE 

GO  TO  70 

C  REVERSE  INTEGRATION 

51  CONTINUE 

00  6C  J=1,N 
JJ  =  JU 

X  (NPT,  JJ»  -  t  (NPT.l) 

C  X<NPT, JJI=O.C 

bc=x(npt, jji 

DO  60  K<=2,NPT 
K=NPT.1-KK 
Kl  =  tf*l 

BD=90*ai*X  <K,  J1 
» (K.JJI-9D 

60  CONTINUE 

IF( IBIAS.  Q.C  >GC  TO  62 
I P  W  R= 1 61  A  -1 
00  61  KK=. ,NPT 
T I M  £=0  T’KK 
K=NPT* 1-K« 

61  X(K.NP21=TIME»»IPWR 

62  CONTINUE 
GO  TO  70 

C  GENERATE  UNIT  DELAYS 

91  CONTINUE 

00  93  1  =  2, NP1 
11=1-1 
X  (1 ,11  =0. 

00  93  K  =  2  ,N°T 
K1=K-1 

93  X<K,I»=Xm,Ill 

GO  TO  «1 
70  CONTINUE 

SC=1.C 

00  80  1=2, NP1 
SC=SC*OT 
DO  8C  K=  1 , NPT 

80  X  IK, I)  =SC*X  (K,I» 

81  CONTINUE 
IFUPR.LT.AIGO  TO  99 

IF( IROUNO.EC. 1»  »RITE(ILT,17F) 

IF< IROUND.EQ. 01 WRITEJ ILT, 179> 

00  180  1=1, NP2 

180  WRITE  1ILT.1 1C  M  X  I K  ,  I  J  ,  K=  1 ,  NPT  > 

C  WRITE<ILT,lt 

99  CONTINUE 


74 


c 

Clio  F0RMAT(M1X,F12.  El) 

110  FORMAT  (20  (IX,  F5. 21  I 

178  FORMAT (10  .InHNOISY  X  MATRIX! 

179  FORMAT C1GX.8HX  MATS! X 1 

I  fcrmat</» 

REFORM 

ENO 

SUBROUTINE  FIXCG, F.C.D.X.N.NC.SIG.NDIP, IEIX I 
C  - 

c 

C  ESTIMATE  NOISE  INTENSITY  SIG  (ASSUME  WHITE  NCISE> 

C  CORRECT  NdSY  rt  AT  R IX  =  C 
C  <P>  DENOTES  NOISE  MATRIX  FO°  UNIT  NOISE 

C  NC  IS  THE  NONZERO  SUBMATRIX  OF  F  =COV  OF  NOISE 

C 

OIMENSION  GCNOIM,  ll.P  (NDIM.  II  ,C  CNDIM.  II  ,0  (NOI*.  II  .  X  (1 1 
IFdFIX.EC.OGO  TC  51 
JCT  =0 
SIG=0.0 

3  SUMCET=u.O 

CALL  GxO0CT(G.0.G0ET,X,0,N,NDIM,0> 

JC  T  =  JCT  *1 

IF ( JCT • EO • 1 1 OET  G=  GQCT 

00  5  J=1,NC 

no  7  11=1. N 

00  7  JJ=1.N 

C(II,JJI=G(II ,JJI 

rr  (JJ.EO.  J!CdI.JJ)  =  PUI,JJ* 

7  CONTINUE 

CALL  GXRDCT(C.D.0r T.X.N.N.NOIM.O* 

SUMCET  =  SU(QET  *OET 
5  CONTINLE 

SI=GDET/SUMCCT 
51  CONTINUE 

no  9  1=1, N 
00  9  J=1,N 

9  C(I,J)=G(!, Jl  -SI*F(I,JI 

IF  ( IFIX.F1.  Cl  GO  TO  11 

CALL  GFROC  T  (C.O.CCET.X,;  .N.NOI’I.C  ) 

IF CCCcT.LT. G. C1SI=SI/2.C 
IF  (COST. LT  .  0.  Gl  CO  TO  51 
IP ( JCT.GF.5IG0  TO  11 

IF  (COET/OETG.  GT  .0  .  1>  CALL  MEQUAT(N  .N.G.C.NCI  M.  II 
SIG=SIG*ST 

IF (COET/DSTC. GT  .0 .1IGC  TC  3 

II  PET  URN 
ENO 

SUBROUTINE  eUILOR < A , X , N, M A  X  I 

C  — - - - - 

C  CONVERSION  MATRIX!  REVERSE  INTEGP  AT IOF  —  I.R.  MODELING 

OIMENSION  A(MAX,1I,X(1I ,Y(2:i 
OOUBUE  PRECISION  ct.y 

COMMON  ZOAO /ISPN,XMSB,0ELTA,SIG2,CT,0I.flIAS.19I AS 
COMMON  /IO/IR.ILT.IPR.IRCUND 
NM1=N-1 
00  11  1=1, N 
vm*b.3 
00  11  J=1,N 
11  A(I,JI=3.C 

A(N,N»=1.C 
C 

00  20  J  J= "  «  NM  1 
J=N-JJ 
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00  15  KK=1,2 
K=KK-1 

00  15  I=J.NM1 

15  A  «I*K.  J1  =  A<  I*  K»  J»  *AII*1.  J*l>*  <1.G-K-K» 

21  CONTINUE 

00=1.0 

00  22  J=2.N 
00=00*01 
00  22  I=J.N 

22  AII,Jt=QO*ACI.J> 

C 

00  25  1=1. N 

x  F  (  IPR,G£.3)»nnc  111  T,  5  J  |A  » I.  JJ>.  JJ=1  »N* 

00  25  J=1.N 

25  V  <1  »*V  (I»  +A  (I  t  Jl**  <J> 

00  28  1=1. N 

25  x<n=T(i»/rci> 

IM  IPR.GE.3)MRITS  <6.71  <X  CI>  ,I=1,N> 

5  FORMAT (2X.1QG12  .51 

7  FORMAT!*  ESTIMATEO  PARAMETER  VEF TOR* • / • 1 CG 13 .6 ( 

RETURN 

ENO 

SUBROUTINE  GKftOCTIX.V.OET.XLAFOA.NN.N.MAY.I  Cl  T) 

C  - 

DIMENSION  X  LA  FO  A  <  1 ) 

DIN  ENS I ON  X<NAX,1».X<HAX,1I 
DOUBLE  PRECISION  A.e.C.O.E 
INTEGER  NLM<2.2d 
OOUELE  FR-CISION  ct.sc.ac.rc 

COMMON  /OAO/ISPN. XFSN.0cLTA.SIG2. CT.OI.BIAS.IOIAS 

COMMON  /IO/IR.ILT.IPR.IRCLNO 

IGKR=1 

IFJN.Nt.HGO  TO  3 

V  Cl. 11=1. 

DET=X(1,1» 

GO  TO  61 
3  CONTINUE 

00  6  1=1. N 
DO  6  J= 1. N 
6  V  (  J,I»=XCJ,I» 

A=  1  ,C 

00  M3  1=1. N 
0=0 .0 
L  =  I 
M=I 

c  FIND  LARGEST  ENTRY  AIL.M)  IN  TMELOHEP  DIAGONAL  SUBM-T-'IX 

c 

00  15  J=I.N 
00  18  K=I.N 

IF < ASS « Y <K, JI I .LE . DABS <91 ) GO  TO  18 

e--Aes(T(K,jn 

L=K 

M=  j 

18  CONTINUE 

C 

C  INTERCHANGE  RONS 

c 

IFa.EO.IIGO  TO  2m 

00  23  J=1,N 

c=y  a.  j» 

Y(L.J»=YII.JI 
23  Y4I,J»=C 

C 
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c  interchange  colunns 

c 

24  IE (M.EQ.II GO  TO  29 

00  28  J=1,N 
C  =  Y  <J,rl 
V(J,MI=Y(J,II 

28  Y(J,TI=C 

29  NUH(i,II=l 
NUH  (2,  I »  =  1 
8=Y  11,11 
Y<I,I>=A 

DO  42  J=1,N 

IFIJ. EQ.IIGO  TO  42 
C=-  Y  (I  ,  J) 

Y(I,JI=j.J 

00  41  K=1,N 
0=Y(K,II'C 
e  =  Y  <K, JI»B*0 

if <ca3s <e i .lt .1.00-1 :»OAes  <oi ie=: .; 

41  Y(K,JI=E/A 

42  CONTINUE 

43  A  =  3 
C 

C  RESTORE  COLUNNS 

C 

DO  58  1=2, N 
J=  N  ♦!- 1 
K=NUM<2, Jl 

IFIK. EO.JIGC  TO  52 
DO  51  L=1,N 

C*Y  (K.LI 
Y (K,L)=Y( J.LI 

51  Y<J.LI=C 

52  K=NUH(1,J1 
C 

C  RESTORE  ROMS 

C 

IFIK.EQ.J1G0  TO  58 
DO  57  L=1.N 
C=Y  (L,  81 
Y<L.KI=V(l, J| 

57  v (L, Jl =C 

58  CONTINUE 
OE  T  =  A 

61  CONTINUE 

IF ( ICPT.NE. II  GO  TO  103  0 
IF  mi, II  .LT.O.CIGC  TO  1C5C 

C  . 

SC=1.C 

00  200  1=2, N 
SC=SC*OT 

IF  <  IGKR.EQ.  01  XLAHOA(II*Y  Cl, II  7Y  (1,11 
IFaGKR.EO.OI  GO  TO  2CJ 
A  =  Y  II,  II 

IF(Y(I,ri.LT.0.0» A=3 . c 
IF <  IGKR.EQ.  2)  A*  ABS  <Y  (1,1  II 
XLANOA (II=OSOPT  (A/Y(l, II  | 

C  IF  (Y  <1,11  .LT.  £,.  0»  XLAMDAIII =-XLAMDA(II 

XL  ANOA  <II  =  SC*  XL  ANOA(I) 

200  CONTINUE 

xlanoa<ii=i.odo 

IF (IPR.GE.l)MRITE (6,1061 (XLANOA(I) , I=i,NI 
NPP  =  N 

IF<I0IAS.NE«OIN OP *N-1 


77 


CALL  BUILOR  IT  .XLAMCA.NPP  ,)>AX> 

136  FORMAT (5X.*SVNTHETIC  PARAMETER  VECTOR*,/, 10 G1 2. 5) 

10CC  CONTINUE 

RETURN 
ENO 

FUNCTION  :OHB  IN , M  ) 

C  CALCULATE'  COMBINATION  M  CUI  OF  » 

IFIN.LE.OIGO  TO  99 
L  *  1 
LO  =  1 

IF(h.£O.C>GO  TO  8 
HN1=N-M*1 
00  5  I *HN1 , N 
5  L  =  L*I 

00  7  1  =  1,  •* 

7  LO=LO*I 

8  CON  8*L /LO 

99  RETURN 

END 

SUBROUTINE  BU  ILOZ  I  Z,», NP1 , NPT  , N Cl M ,  (.FIX) 

c  - - - - . . 

D I  MENS  ION  Z  INDIM.l), Rll) 

OOueLt  RftrCISI0N  CT 

COMMON  /D'  O/ISPN, XMSB,0ELTA,S1G2,CT,QI,9I AS  ,191  AS 
COMMON  /IO/IR ,ILT, IPR, IRCLNC 
TIME=0T*N°T 
I Cp  T  =  NFI X  *1 

GO  T0(2D1,1S1«1G1«2C1)*XCFT 
101  SC=  OT 

00  24  1=1, NP1 
ZII-,  1*=SC»NPT 

IF(I.GE.2)Z(I,2I=0T*SC*C0MB(NPT-1*I,I) 

24  SC=SC*OT 

00  166  J=  3 , NP  1 
00  166  I  =  J,  NP  1 
A8=TIME»* II+J-?) 

AC=1.0/(I*J-3I  -  l.C/II+J-2) 

166  Z(l,JI=A8*AC»0T 

Z (3,3>=Z(3,2) *0T» (NPT*1 .01/2. C 
WRITE (ILT  ,16C > 

GO  TO  331 
231  CONTINUE 

00  213  J=  1 ,  NP  1 
00  213  1  =  J, Np  1 
IF(I.EQ.1»Z<1,H  =  TIME 
21 C  Z  (I,  J)  =  (TIME*  *11+  J-l  >>/(!•*  J-l) 

WRITEIILT.iei ) 

331  CONTINUE 

00  174  1=1, Npl 
DO  168  J= I , NP 1 
168  Z(I,J)=Z(J,H 

IF (IPP.GE.3)HRITt  <ILT,22G)  I Z 1 1 , J) , J=1 , NP1 ) 

174  CONTINUE 

16 C  F  CRMAT ( 10  X, *  OUANT,  NOISE  *> 

161  FORMAT  (1CX,*0IAS  effect*) 

220  FCRMAT<?X,5(2X,G13.6M 

return 

END 

SUBROUTINE  PPTHAT I  A , M , N, NO IM, LOCI 

c  - - 

c 

C  PRINTS  A  MATRIX,  ANC  AN  IMEGER  (PERHAPS  A  LOCATION)  IF  LCC.S 
C 

DIMENSION  A(NDIM,1> 
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IF!LCC.GE.1IWRITE!6,5»L0C 
DO  31  1=1, N 

31  WRIT£!6,i*MA  (T,J»  ,J=l,i4> 

5  FOR«AT«*  LOCATION/ INTEGER®* ,I5> 

1A  FCRMATI2X.10G13.6) 

RETURN 

ENO 

SUBROUTINE  RESPON ( X, V, N, GAMMA ,XLAMDA,MPil 

DIMENSION  X!U,  V(1I.GAMMA!1»,XLAMCA<11 

DOUBLE  PRECISION  XSAV.AO.BD 

NM1=N-1 

NP1=N+1 

NPNP1=N*N*1 

NPNF2=N*N*2 

OC  IB  I=1,NPNP1 

19  XL AMDA 1 1(  =9  .0 

xsav=o 

00  2C  K=1,MP1 
IFIN.E0.1IGC  TO  25 
00  21  1=1, NM1 
J=NF1-I 

21  XLAMOA(JJ=XLAMOA!J-H 

25  CONTINUE 

00  22  1=1, N 
J=NFNP2-I 

22  XLAMCA!J)=XLAM0AIJ-H 
XLAMDA(1)=XSAV 
XLAMDA  f NP  X I =tf 1*0 
XSAV=a .0 

00  23  I=1.NPNP1 

23  XSAV=XSAV-GAMMA (I*1)»XlAMCA  (I > 

IF ( DABS ( XSA V)  .GE.1.0E10I XSAV=C,3 

20  X(K»=XSAV 
RETURN 
ENO 

SUBROUTINE  VE  OU AT  <NPT , Y, X, NPUL ,  IOFT) 

C  - - 

C  IOPT=a  SET  Y  TO  7EP0 

C  1  OR  2  SET  Y=  X  tPRINT  IF  2) 

C  3  S^T  Y  =  Y  *  CONI T  X 

C  9  SET  Y  TO  ZERO 

C  10  STT  Y=IMPULSE 

C  11  SET  V=STEP 

0IMEN3I0N  X  (1 l,Y(ll 
IF  <  IGPT.S1.  Cl  I0FT  =  9 
00  33  K=1,NPT 

IF ( I0PT.E0.1.0».ICPT.EQ.2IY(K»=X!tO 
IF! I0PT.E0.3>  Yl X»  =  Y(K»*X!i* 

IF!I0PT.G-,9! YIKJ=O.C 
IF  !  IOPT.EO.il  •  Y  IK  )  =1. 0 
33  CONTINUE 

IFU0PT.E3.2I  WRITE  16,6) I Y I K> , X= l, N°T» 

6  FORMAT  !2X  .10G12.5) 

IFIIOPT.EO,  10 )YIU  =1.0 
RETURN 

ENO 

SUBROUTINE  MEOUAT!M,N,B, A ,NDI M, 10  FT) 

C  . . . 

C  IOP  T=0  SET  9  TO  ZERO 

C  IB  EQUAL  TO  A 

C  10  B  TC  IDENTITY 

DIMENSION  AINDIM,1I,BINCIF,1> 

00  33  1*1, N 
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o  O  C*>  <D  oooooo 


00  33  J=i,N 

IF(IOPT.NE.l)  8(1, J)=0,C 
IF(I0PT.£0.1C .AND.I.EQ.J)B(I, j)=i.: 

IF(IOPT.EQ.l)  B(I, J»  =  A(I, J) 

33  CONTINUE 
RETURN 
END 

SUBROUTINE  PLOP(NPT,nf,Y,NDIM,T:,CT,lABEL .INr£P, 


NPT =NUM3  OF  TIME  FTS  (NARNINGl  NCIl  SHOULD  «-.C' 
NFsNUMBER  OF  FUNS 

Y(K,  )  DA"  A  ARRAY  CF  DIMENSION  NCIr'..NF 
TO=INITIAL  TINE,  OT=TIMS  INCPtMENT 
LABEL,  INOEP  =  TITLES  FOR  Y  AND  X  AXES 
DIMENSION  Y (N  CIM, NF| , YYC2) .LAPEL (1) , INOEP  (1 1 
DIMENSION  X (5 12), I8UF(S12> 

COMMON  /IO/IR.ILT. IPR, IP CLNO, TPLT 
M=NF»NCIM 
M 1  s  M*  1 
M2=M*2 
NPT1=NPT*1 
NPT  2-NPT ♦ 2 
X ( 1 1 =T  0 
DO  9  <C  =  2,NPT 
X(K)  =  X(K-1)«T 
DO  8  1=1, NF 
no  8  K=NP71,N0Iv 
Y(«,I)=Y(NPT,  I) 

INITIALISE  mo.  INF,  12In.FAPSF>  MA x. L  E  Nr,7N=e  : : 

call  PLOTMX  (o  0  .  0 ) 

C  SET  ORIGIN 

CALL  PLOT(C.,-. 5,3) 

CALL  fact  R  (5  ,0/6.5) 

C  BEGIN  PLOTTING 

CALL  SCALE (X, 6. 5, NPT, l) 

call  scals<y«i,i) ,10. :,m,i> 
call  AXIS<C..  C.  .UPTIME  (SEC.), 

•-16.6.5, j.,X(NPTl) ,XtNPT2)) 
call  AXIS'O., G. ,16H»ESP0nSES  Y  ,y, 

*16.1C.,9S., YIM1),v<m2» j 
HRITE(6,6)X(NPT1) , X ( NPT2 ) 

HP  I  TE  (6,7)  Y  (*  1)  *  Y  (M2) 

6  FORMAT  (IX, *TG  ,01V  (6.5  0  I V )  *,  A  ( IX  ,  F7 . 3)  ) 

7  fqpmAT (IX, »Yt  ,01V  (10  DIV)*,4  (1X.F7.3)) 

00  10  1  =  1, NF 

Y(NFT1,I)  =  Y(M1) 

Y (NPT2,I)=Y (M2) 

IF(  I.  EC.  1.  OR.  IPLT.EO.  3)  CALL  L  INE  (X,  Y  (1,  I)  ,N  FT  ,1 , 1 
IF ( I.E0.2. ANO.I FLT.EQ. I) CALL  CASHLN(X,Y (1 ,2), NPT, 
10  CONTINUE 

CALL  PLOT (1C. ,3 ..-3) 

RETURN 

ENO 
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USPEC 


“US °EC“  P0-.  R  ADr 

EV4lUATES  THE  MAGNaTUOE  S»rCT«UM 

of  THE  FOLLOWING  PULSE  INPUTS  (SEE  SECTION  4 

OF  THE  REBOOT  I  I 

1,  TR*,-(TI  ONE  CYCLE  TSI  WA VE,  FCLLOWEO  ?v  ZERC  LEV'LL 

2.  TPMT)  ♦  VE  w  (:  L r  CYCLE  TPI  WVA£,  FOLLOW-::  I Y  Z?r  0  LEv'L 

3  .  CEXJ(T)  CSOTLLjTCVY  FULSE  (SO  DECAY),  F  CLLJ  wEU  BY  7E\C  LEv/El 
CEYKTI  S~mE  is  T,  with  =  >=>.  CiCAVTKS  ENV.  -  ONE  TIB:  C0'.ST 

5.  OEXKTI  SAW-;  as  3,  WITH  E*P.  DECAYING  Env.  -  T«C  TIME  C0MS7 

6.  OTP  ( T I  OSCILLATOR  FULSE  WITH  T#I  EN V.  ,FCL LO W 0  BY  Z  EEC  LEVEL 

FO-  PULSES  3  TO  6  The  OSCILLATION  IS  A  COSINE  WAVE 


DIMENSION  Nil  C>  ,F  <2u»  ,FR  <2C  >,  >(20,101 

DATA  N(H/2C/,N(2l/i?C/»M3l/20Z/»NT4»/lCGC/.ME»/2;i:/ 
OATA  F  (l  >/.  Cl  /  .F(  2)/  .  1X,F  <3)/  .5/,  F  .  9C  oi/,  FJ  5)  /  .  990;  9,/ 

OATAF(pl/.99M/,c(T(/i.;i/,F(/>l/:,l/,F(,)/l,j./,  FCi;  l/i.  „•/ 
OATA  F(llt/3. »/,F( 12 >/i:./,F( 131/11./ 

PI=,.Q*AT'  N  (1.0) 

012=2. 3*01 
NN  =  5 
NF=13 
IS1=1 
I S  2  =  6 
N*l  =  3 

00  3  C  J=1,NF 

IF (F(J> ,LE. 1. C) NF 1=NF1*1 

30  CONTINUE 

I TEST=C 

IFCTES'  .NE  .1)00  TO  29 
n (ii =  10: 

F  (1)  =0.1 
I  Sl=  3 
I S2= IS  1 
NN=1 
NFF=1 
100  =  3 

29  CONTINUE 

IF(ITEST.NE.l )IPP=0 
00  lOu  ISTG=I SI « IS  2 
NFF=NF 

IF (ISIG.LE. 2) NFF=NF-NF1 
00  95  1=1, NN 
00  94  J= 1 , NFF 
F0{  JI  =  F<JI-l.f. 

IF (ISIG.ET.6. ANO. F (Jl .LE .C. Ori) FR ( J)=-C ,0 
IF  (  ISIG.GE.  3.  OR.  ITEST.EQ.DOC  TO  31 
FP(  JMF(J»NF1)-1.0 

31  CONTINUE 

GO  TO  (41.41, 51,51,51,61, 71), 1SIG 
41  CONTINUE 

9Fr(0,125r'-2l*N(I) 
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IF<IPP.GE.3>wFIT::lfe,e  >3F 
8F=8F*FS(J1 

IF (IPP. G£. 3 > Tf (6,61 *F 
TWETA=PI*PF 

IF  (  ISIG.EG.2)  Tr4£TA  =  2.0*Tt-ETA 
IF (  Ip- .GE.  2  (WRITE  (6,7 > ThETA 
xx =1.' 

IF  (  THETA.i.E.C  I  XX=STN  <TWE  TA»  /TH£TA 
IF ( IPP.GE.2IWPIT-  <6. 81  XX 

X  y=**»*< 

IF  «ICx.r,E.3»uKITf  (6,  9»  XX 
IF<ISIG.EI.lt  0=SIN (2. C’TPET A) 

IF  «  ISIG.FQ.lt  XXi:y  X*  A  3S»0» 

GO  TO  93 

continue 

A=G  .0 

8=o.o:5»n  m 

IF  (ISIG.EC.M  A=2C3.C/MI> 

IF  ( ISIG.:0.5>  A=<.Ct.S/N«I> 
EA9=£XO(-fc*0) 

IF  (IPP.GE.21WPIT--  (6,l3t-.6,EAe 
XF=0.E 

XI»5.C 

OC  57  XF=i,2 

FF  =  FR  (  Jt  ♦  2.  C*  (KF-ll 

Wip  I?»FF 

C=  CCS  <  W*9t 
S=SlN(W*B) 

IF(IP-.GE.2>WFITE(6,l»>h,C,S 

VR*1.G-E43*C 

Yl=-EAB*S 

ZP=A»YR.W»YI 

7I=A»YI-H*YP 

Q=A»A»-W*W 

yP=yP*Z3/. 

x:  =  xnzi/o 

CONTINUE 

XX=Sa-’riX^«XF  +  XI«yII 

IF ( IPP.SE, 2tN  FITE  lt,:6>XF,XI, XX 

IF< j.EQ.nfiixc=xx 

go  TC  33 

CONTINUE 

9=0.0C5*N(I> 

IF ( IPP.GE.2 )H  FITE  «fc.?3l e 

xp=:.; 

XI=3.2 

00  67  XF=:,2 

FFsFCi Jt+2.c* (<F-1) 

W=PI2*FF 
JmETA=pi»ff*R 
0*1.  C 

IF (THETA.NE .0 . C 10=SIN (THETA l/THETA 
C=C0S(T9£TA I 
S=SIN( TmETAI 

IF  ( IP7.GE.2 IWfcITE (6,2-1 W.TMFT A,D. C.S 

YP=1.0-O*C 

YI=C»S 

XP=XA.YI/N 

XI  =  XI- VR/W 
CONTINUE 

XX*SOPT(X»*XP*XI*XI» 

IF  (IPP.GE. 2 INRITE  (6. 15 MR, XI,  VX 
IF (J,cO.Ncl)yC*Xx 
GO  TO  93 


roNTiuut 

C3NTINU; 

X TJ,I»=XX 
CONTINUE 
DO  9c  J=1,NFF 
XX=X(JtI> 

TF(ISIG„G".3. AND.Mrr,Gi.NFl)X>=x<J.I»/x: 

x  ( j.lt  =-2:0 

IF (XX.Gr .1.5c -1C » X (J,I >=2G. C» AtOGlu (XX* 

COWTIMU- 

CONTIUUi 

WRITf.  C6.11IISIG,  <N  (I  t  ,  1=1,  N‘N> 

WRITt(6,i» 

00  9?  1=1. NFF 
FRO=FP  <I( 

IF ( ISIG.G'.3> F9Q=cc( I)+1 ,c 
WRIT£<6.9(FSQ,  «X<I,J>  ,J=1,N*;> 

*RITE<6.1> 

CON  TlNUi 


FORMAT  </l 

FCRMA7f?X.F12.3,9Fl2.1) 

F0KP4T<2x,*THhT4=»,Fi3.e* 

F0O«4T (2X,*Y=SlN(  »/(  >=*,F1T.6> 

FORMAT (2X,*V. Y=*,F13.c> 

FORMAT  <10  *,*m  AG.  IN  rl°  FOR  01  FFEn,  FNT  FRcC  «N  FT  »,  *  ISI&“ 

♦  Il.//,7X.*F/Ft/»,1£(7X,I4,,1XM 
FORMAT  <14  x,  |.'|  I*  ,  NX  »  » 

FORMAT<2X,*A=*,Fi;.5,‘>  c  =  *»F1C,5»*  t  A8=*,  Fli..  5» 
FCRMAT(?X,«M=«,fi:.5.*  CCS=-»,Flf  .5,*  SI  N=  »,  Fj  £  .  i  > 
FCRMAT<2X.»  X*  Rf Al=»,Fi: .  =  ,*  »*  I  f  A  G=  *» iF 10 .9  I 

FCRMAT<2X,*e=*,Fl£.F) 

FORMAT  (2X,«W=  *,fi  :  .5,  «  TmFTA  =  »,fi:  ,5,»  SIM  >/•  >*.F 

♦  /,2X,*COS=*,Fl3.‘>,*  SIN  =  *,FU.5) 

STOP 

FNO 


*'-&r 
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5  MISSION  b 

5  °f  b 

§  Rome  Air  Development  Center  | 

n  RAPC  plans  and  execute*  research,  development,  test  and  v 

»  selected  acquisition  programs  In  support  of  Command,  Cont/Lol  a 
«  Communications  and  Intelligence  (C3 1)  activities.  Technical  ^ 
%  and  engineering  support  within  area*  of  technical  competence  Q 
fc  is  provided  to  ESP  Program  Office*  ( PO* )  and  other  ESP  § 

D  element*.  The  principal  technical  mission  area*  are  S 

*  communication*,  electromagnetic  guidance  and  control,  *ur-  a 
R  veitlance  of  ground  and  aerospace  objects,  intelligence  data  v 

*  collection  and  handling,  information  system  technology,  v 

|  ionospheric  propagation,  solid  state  sciences,  microwave  2 
Ig  physic*  and  electronic  reliability,  maintainability  and  Q 

h  compatibility.  <& 


